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A  Locally  Linearized  Panel  Method 
■for  Calculation  of  Tr  i -dimens  i  onal 
Wing  Oscillation  Pressure  Distribution 
at  Trans/subsonic  Speeds 


Tong  Binggang  Zhuang  Lixian 
<Dept.  of  Modern  Mech.,  Un  i  v .  of  Sc  i  .  6t  Tech  of  China) 

Li  X i an  1  i  n 

(Computational  Center,  5th  Inst,  ot  Ministry  o-f  Aero.  Studies) 

ABSTRACT 

In  this  paper,  the  theory  of  disturbance  potential  -flow  is  used  as  the 
basis  for  a  local  linear  hypothesis  to  construct  a  tr i -dimensi onal  down-wash 
integral  equation  for  wing  oscillation  at  trans/subsonic  speeds.  The  method 
for  calculation  of  the  trans/subson i c  flow  Kernel 

function  is  expounded  upon.  The  Generalized  Doublet  Lattice  Method  is  usea 
to  calculate  unstable  atmospheric  pressure  distribution  over  the  rectangular 
resonant  cavity  of  an  aircraft  wing.  The  results  indicated  the  local  method 
compared  accurately  in  calculating  nonlinear  characteristics  of  unstable  flow 
forces  at  trans/subsonic  speeds.  In  regards  to  calculation  of  flow  forces  at 
supersonic  speeds,  more  research  is  required  to  find  a  method  of  numeric 
derivation.  1  J 


I .  Foreword 

* 

As  everyone  knows,  potential  linear  flow  is  often  described  by  means  of 
a  partial  differential  equation  of  second  order.  If  the  equation  is  linear, 
the  fundamental  solution  is  found  using  repetitive  addition.  Boundary 
conditions  for  the  physical  plane  may  be  expressed  in  the  form  of  an  integral 
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In  the  formula,  S  is  the  solid  boundary  of  the  field  of  flow,  added  to  the 
interruptive  residual  surface  of  the  tail,  V „  is  the  velocity  cerivec  at 
plane  S.  K.  is  the  fundamental  solution  to  the  partial  differential  equation 
of  tne  i  degree.  dS  egresses  a  possible  fundamental  solution  for  pane! 
weight  of  dS.  and  V n  are  pre-de term i ned  functions  derived  from  the 

formula  and  surface  Planar  conditions,  is  a  possible  position  within  the 

selection  order.  The  remaining  elements  deal  with  certain  known  functions. 
After  resolving  C..  distribution  functions,  the  field  of  velocity  and  fieid  of 
pressure  can  be  expressed  integrally  as: 


D,K.dS 


which  gives  us  a  true  derivative.  For  the  spec i a  :as ?  of  conditions  of  the 
lift  plane,  we  can  also  obtain  an  integral  equation  for  extracting  load 
di str i bu  t i ons . 

The  Panel  Method  of  calculation  is  one  o*  the  most  important  in 

[1] 

hydromechanics  for  finding  numeric  values.  The  af oremen t i oned  integral 
equation  can  be  used  to  extract  an  accurate  representation  of  flew  since  it 
traverses  opposing  boundary  conditions,  curved  surface  ooundaries.  and  ode 
distributions  wnich  begin  to  depart  from  theory.  Because  the  definition  of 
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the  realm  of  the  integral  equation  has  already  degenerated  to  that  of  a 
curved  surface,  tne  computations  derived  using  the  integral  equation  as 
compared  to  direct  extraction  through  a  differential  equation  < Cauchy's  Test 
for  example)  fell  short.  As  a  result  of  these  problems  established  in 
integral  equation  (1),  the  Panel  Method  proved  far  superior. 


The  Panel  Method's  suitability  for  solving  problems  of  potential  linear- 
flow  under  normal  and  abnormal  circumstances,  both  for  non-compressed  and 
compressed  flow,  has  given  rise  to  its  extensive  and  successful  application 
C2,  3,  4,  53.  Speaking  in  regards  to  small  perturbation  at  transonic  flow, 
although  it  has  already  been  proven  that  the  potential  flow  theory  is  a  very 
good  approx imat i on ,  however,  because  velocity  flow  equations  are  nonlinear  in 
nature,  direct  application  of  the  Panel  Method  proved  to  be  very  difficult. 
This  paper  used  the  various  proposals  set  forth  by  Spreiter  on  a  local  linear 
theory  16,  7,  33  to  construct  a  panel  method  for  applications  of  unstable 
transonic  flow  which  is  called  the  "Locally  Linearized  Panel  Method."  In  the 
selection  of  Lissajous's  Dispersion  Model,  this  document,  in  regards  to 
trans/subsonic  flow,  proposes  the  use  of  a  Generalized  Doublet  Lattice 
Method.  This  metnod  was  used  for  all  experimental  calculations  in  order  to 
clearly  demonstrate  its  effectiveness.  As  for  the  criticality  of  transonic 
flow  <to  include  the  potential  flow  in  the  presence  of  shock  waves)  it  seems 
rational  to  use  Lissajous's  Figures  as  a  model  for  construction,  but  further 
studies  are  required. 


II.  A  Down-Wash  Integral  Equation 
for  Wing  Oscillation  at  Trans/subsonic  Speed 


-3- 
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In  recent  years,  attention  has  been  towards  estaol  ishment  of  a  method  *or 
a  Kernel  function  at  tr  ans/ 'sub-son  i  c  speeds,  in  the  western  world  by 
Cunningham,  A .  M .  [?,  IS]  am  _.u,  D.D.  £11,  12,  13]  et  al  .  Cunningham  s-rr,:'. 
does  not  use  an  integral  equation  based  on  linear  theory,  but  from  within  the 
Kernel  function  he  makes  several  adjustments  in  regards  to  Mach  number 
and  frequency  K.  _iu*s  method  leads  into  transonic  speed  ac  c :  '  e  •  a  .  :r 
force,  and  attempts  to  establish  a  velocity  acceleration  force  equation.  The 
next  step  leads  to  a  transonic  speed  down-wash  integral 
equation.  However,  we  c i soever  so  n ; s  a.:,  at  ions  were  in  erre" . 

Resultingly,  a  solution  remains  unestablished. 


As  a  starting  point,  this  paper  devises  a  local  linearized  hypothesis 
from  the  force  velocity  equation  for  unsteady  state  small  perturbation  and 
derives  a  down-wash  integral  equation  for  wing  oscillation  at  transonic 
speeds.  Consider  an  arbi trary  planar  surface  of  non-torsi onal  quality-with  a 
given  thickness  as  the  wing.  By  placing  it  in  the  surroundings  of  a  near 
sonic  horizontal  flow  influx,  the  horizontally  positioned  planar  surface 
serves  as  a  small  scale  resonance  cavity,  as  the  vibrations  reach  near 
stable  conditions,  the  small  perturbation  force  velocity  of  flow  can  be 
expressed  as: 


3> 


In  this  equation,  +  expresses  the  wir.g  thickness  that  gives  rise  to  steady 
perturbation  force  velocity,  r-:  •  expresses  the  wing  thickness  that  gives 
rise  to  unsteady  perturbation  force  v e 1  pc > t y . 


Based  on  the  smal 1 


! 
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perturbation  -force  velocity  hypothesis,  we  have  *'  •  o  ,  ,  arc  tr*  s  allows 

I 

I  us  to  -formulate  the  approximation  scale  equation  be  lew: 


( 1  —  A/  i )  o„  -6,,—  <&,.  =  Q 


M'a1  T 

<P..  Q  [11 


x~L~u — 


Here,  U  and  M  differentiate  the  inward  velocity  and  Mach  number  that  is  the 
most  recent  or  ixmediate  Mach  number  value. 

( 

i 

By  takinQ  the  thickness  distribution  of  the  wino  and  the  vibration 
displacement,  the  differentiation  can  be  expressed  as: 

z,-h(x,y)  ft  i.  =  /(x,y)3'"  <7> 

and  with  equation  (4>  boundary  conditions  are  expressed  as: 

( *L\  .c/* 

•  ox  <s> 


while  equation  (5)  expresses  the  boundary  conditions  as: 
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Herein.  k-aLlU  u  the  conversion  frequency.  L  is  the  reference  length. 
Equation  (.4)  is  a  nonlinear  eauat  i  on ,  equation  (5>  is  a  variable  coefficien: 
in  linear  form.  In  order  to  derive  these  equations,  we  used  the  locally 
linearised  hypothesis  in  all  instances.  The  problem  of  boundary  limits  in  a 
differential  equation  became  a  related  problem  in  integral  e : „ a  c i on  s .  so  t  r  e 
Panel  Method  for  derivation  of  numeric  values  was  once  again  utilised.  In 
regards  to  the  local  linearized  Panel  Method  of  equation  <4>,  you  may 
reference  footnote  C14]  which  discusses  in  detail  a  Panel  Method  of 
calculation  for  unsteady  state  trans/subsonic  flow  oscillation. 


According  to  the  locally  linearised  hypothesis,  we  can  take  tne 
coefficients  from  the  various  known  functions  of  equation  <5>.  Afterwards, 
at  the  proper  time,  we  take  these  parameters  and  once  again  use  cnem  a= 
values  for  the  originally  known  functions.  We  first  discuss  A/,<1 
pre-cond i t i ons  for  establishment  of  the  trans/subson  i  c  speed  flew.  IJ  :h 
Goe  t  n er  t '  s  t  r  an  sf  orrr.a  t  i  cn 


x'~x,  y'^0,y,  z’-0,n  <p~<ple•’ 


wherein 


0i»\/  I  -M\  , 


_ u»A/‘  ,  A 

a~~u jr'rtT’ 


t_M' {•/->- 1)6„ 

A  u 


•'  i : 
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and  both  become  constants.  Hence,  equation  <  5  >  becomes 


In  this  manner  , 


V'» 

i 


dx'*  dy'‘  dy*  ’ 


.  1 


the  hypothesis  -for  p  ,  in  all  respects  o-f  continuity,  can  be  differentiated 
(without  shock  waves)  within  the  Green  formula 


Groer  u«<*.*;fV;, 


r' -  iy' -y\)--r  (y-Zl)> 


*nb  util  i  cat  ion  of  p,  for  g\  «o  Diane  surface's  counter-character  ;  s*  i  cs  ,  we 
Qe  t 


i-jj  &pt<x[,y[)  -JL(-±Jl,,LL'\dx'dy[ 


z ' 31 0  planar  projective  surface  of  the  wing  and  tail  flow  surfaces.  A<P> 
is  S'  second  reading  p,  difference.  Transformed  and  intre pretec  v • a 
e  q  u  a  t : : n  ■  13'.  it  is 


—  7  — 
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<p(x,y,z)  =  -  J-et,'\\±<p(xl,yl'>e',’'~^-——\dx,dyl 


where 


r  =v/  (x-Xi)1- ^‘,Z(y-yi)2~  (z-zj‘2 


0  and  once  again  the  linearized  pressure  expression 


P - pJJe'  a  -L.  ( <pe  o  y 


is  utilized.  We  derive 


*i 


i n  the  form 


Ap  =  p,.»-  -p»->* 


and  equation  <16>  is  substituted  in  equation  (14)  where  we  have 

<?kx,-j,z)  =  -  --1-  ^1) 

■l-TP-C  Jj  vr  \  r  / 

e  . 

.  .  — 4  .1  1  •*  1 

j  ±p</.,y,)e  -■  J/.je’  a  dxxdyx 


(16,' 


(  1  7> 


rakino  the  derivative  z  from  the  above  equations,  we  can  obtain  a  limit  for 
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z  —  -) 9 .  After  performing  the  transformation  calculations,  it  takes  the  form 
of  a  trans/subson i c  down-wash  integral  equation: 


w  =  -  y,)KLt~x.,  y-'ji  <t>„  6,.)dx..  dy> 


Here,  the  Kernal  function  K  takes  the  specific  form  of 


a 


K  (.x,  (/>  =  lim  e~  “  -g~T  j 

a-»t)  W  -m 


(  1 '?  > 


•and  t  =^¥-0)  (y:-z:)  . 

It  is  important  to  point  out  that  the  parameters  a .  a,  k,,0t  ,  etc.  within  K 

have  already  used  the  known  functions  £,y,  .  or.ginally  derived  from  eauation 

(4)  (taking  node  coordinates  x.y  as  parameters'.  They  reflect  that  portion 
of  the  coupling  coefficient  of  trans/subsonic  force  velocity's  steady 
state  opposing  the  unsteady  state. 

In  r-e;..-:s  to  the  problem  of  unsteady  states  at  trans/supersonic  speeos* 
it  is  possible  to  apply  an  integrated  method  for  transformation: 

x  =x,  y‘  =  '0i  V,  01  US') 

and  derive  an  identical  form  down-wash  integral  equation  (IS).  .  The 
difference  >s  in  the  integral  equation's  dependent  doma i n  of  nodal  regions. 
Tne  Kernel  function's  basic  form  must  also  be  altered  to 
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Deriving  the  numeric  values  for  the  Kernel  -function  of  the  integral 
equation  in  the  Panel  Methoa  is  a  crucial  step.  Obviously,  tne  Kernel 
function's  accuracy  in  calculation  <1?>  is  extremely  difficult 
and  we  can  only  make  an  approximation. 

It  is  for  this  reason,  when  sol\ir.g  normal 

problems,  that  we  can  only  figure  as  positions  along  the  aircraft 

wing/s  surface.  Moreover,  we  cannot  calculate  numeric  data  tor  the  tail 
surface  flow.  However,  with  the  use  of  equation  <19>  we  can  derive  the 
necessary  figures  for  normal  tail  flow  as  a  function  value.  In  tne  integral 
equation  ( -oc,x  —  x:)  ere  these  values  u,  expresses  the  x-value  for  tne 
rear  section;.  However .  : f  we  are  speaking  solely  in  terms  of  the  norma’ 
tail  flow  in  correspondence  to  a  particular  point  on  the  wire's  sun-?:*,  the 
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rate  is  i naeeo 


:re  region  of  (-ao,x-x.) 


:  -  *  inteoral 


'  ;  .V/,  to  substitute  for  Mi  .  -ence.  t.-.e  successful  use  of 


the  Kernel  fun  u  : .  zr.  in  tne  integral  equation  at  this  stage  couio  oe  acnievec 


in  linear  theory  calculation  methods  (see  footnote  [5]). 


In  the  Doublet  Lattice  Method,  K  is  normally  altered  to  become 


K(x-xt, 


and  herein  we  have 


k-'.  ir’  1  -T-H»  K'+tt)  ‘-“'l" 
1  -» 


rv—V(y-y  + 


and  is  like 


/  =  /,-/,=  \  (-)Ji-  j  <•••)«/« 


£>  a.  _r  o/i 

*  ’  aTT’  **  ’  Jr, 


where  we  obtain 


£.=  -(  _JL^‘  du+MlJ^L  (r,*0) 

“ '  J  (1-t-w1)  '  z?v/  l-r«; 


so  we  have 


(23) 


(24) 
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i 

f 


(x  -  x.)  —  MR 


01  r, 


(.x-x.Y~8lr\ 


■ 


and  the  numeric  values  o-f  equation  (24)  are  calculated  in  a  very  convenient 


i  manner  [5], 

i 


Now  we  can  use  a  problem  relating  to  stable  characteristics  with  a  low 
•frequency  development  <K  <<  1)  and  take  these  small  small  values  to  use  as 
exponents  -for  our  -function  values.  This  will  give  us  a  start  and  make  the 
calculation  of  equation  (25)  much  more  convenient.  Just  as  we  get  to  the 
third  stage,  we  get 


KtmF(x,)-F{x,),  (r,#0) 


(26) 


where  the  functic 


s  de  fined  as 


and  in  it  is 


Fa)  (*-;.) +rU)r‘ 

&-)] 

+  -r=7 - y.,  -rla[(jt  — A)  — r(A)]l 

2[(x  -a) -rU)_r(/.)  u  J 

«“(x  —  —  8*  1 


(27) 


•zZ  as  r\*Q  .  a  separation  of  the  odd  nodal  points  at  the  rear 
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comparison,  the  charts  found  in  footnote  C15]  are  provided  as  a  contrast  tc 
the  experimental  numeric  data  and  the  results  of  a  linear  theory  approach . 

Figure  1  and  Figure  2  represent  unstable  pressure  curves  wh  e  r  M*<r  . r .  .1 

is  clear-  in  Figure  1  that  when  M=0 . 7  this  method's  calculation  of  ,  C  in 
curvilinear  measurements  and  numeric  data  is  much  better,  especially  along 
the  forward  naif  of  the  wing.  Theory  and  experimental  results  were  in  acccrc 
wit*  ?a:r  other.  However,  when  using  the  linear  theory  to  obtain 
results,  there  was  a  great  disparity  that  was  quite  obvious.  Figure  2 
expresses  a  phase  curve  that  clearly  shows  a  discrepancy  between  theory  and 
experimental  resul ts,espec i al 1 y  for  the  vicinity  of  the  exterior  strain 
measurements  for  3/4  of  the  wing.  Footnote  [15]  makes  a  further  comparison 
and  in-depth  explanation.  It  must  be  pointed  out  that  these  discrepancies 
were  at  M*0.7  conditions,  which  includes  wing  oscillation  and  the  effects  of 
wind-caused  vibration. 

Figure  3  illustrates  I  C7#  J  — .as  a  curvilinear  surface  when  M=0  .  S 

(near  Mach  number  boundary)  fo*  the  local  met nod.  At  the  same  time,  it  also 
provides  conditions  for  M=0.7  and  from  the  experimental  results  o* 

footnote  C151.  For-  both  of  these  represented  conditions,  we  used  averaged 
numeric  values.  Of  course,  this  comparison  is  not  one  of  theories,  but  is 
provided  just  to  serve  as  a  reference  point  because  the  flow  at  M=0.?  is 
clearly  a  shock  wave.  However,  as  a  reference  it  is  very  obvious  to  see  the 
experimental  results  coincide  with  theoretical  for  the 

forward  portion  of  the  wing.  For  the  rear  section  of  the  wing,  newever , 
there  is  a  large  discrepancy  between  theory  and  experimental  results.  This 
is  because  under  these  conditions  is  influenced  by  the  shock  wave 
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of  the  wing.  We  used  “averages"  here  in  order  to  make  a  relevant  comoer 


Figure  1  Figure  2 

Key:  a  -  this  method;  b  - 
linear  theory;  c  -  test  results 

In  summation,  we  can  make  the  -following  conclusions:  this  caper 
proposed  a  local  linearized  trans/subsonic  speed  oscillation  down-wash 
integral  equation  as  well  as  the  Generalized  Doublet  Lattice  Method  for 
trans/subsonic  speed  flow.  It  is  possible  to  calculate  at  small  widths  th? 
small  perturbation  amp  1 i tudes. f or  the  wing's  planar  surface  pressure 
distribution.  It  is  aiso  possible  for  reactions  to  occur  at  transonic  flow 
in  an  indiscriminate  manner  for  nonlinear  forms.  This  is  why  the  method  is 
superior  to  that  of  linear  theory.  If  we  should  use  a  blend  of  <1?>,  (20), 
and  (21)  for  our  Kernel  function,  the  original  data  values  can  provide  a 
numeric  foundation  for  calculation  of  trans/subsonic  flow.  However,  if  the 
theory  could  be  designed  in  a  manner  to  effectively  incorporate  Lissa,ious  s 
Figures,  it  would  advance  us  a  step  further  in  more  profound  research. 
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Figure  3 

Key:  a  -  "reference  comparison"  data;  b  -  shock  wave 
field;  c  -  test  results  for  this  method 

It'  -shock  waves  are  manifest  from  the  wing's  planar  surface,  it  is  necessary 
and  proper  to  include  the  relationship  of  these  waves,  as  well  as  their 
interval.  The  unsteady  state  down-wash  integral  equation  handles  this 
nicely.  In  normal  problem  solving,  it  is  difficult  to  ascertain  under  what 
circumstances  it  is  best  to  use  the  Panel  Method  to  accurate  1 y- de term i ne  the 
point  from  which  shock  waves  manifest,  or  its  type  and  characteristics.  hc 
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Abstract 

u  this  paper,  the  do  *  n-w  asa  in  tec rai  equation  lor  the  smail 
7*:-urbnt:on  transonic  potential  flow  past  an  oscillating  wiaj  is 
established  0£J  * jje  basis  of  local  linearization  concept.  The  method 
c.a : c u i a t : n?  the  tran-/ subsonic  kernel  function  is  discussed  in 
'•  ii  To  Cw.i  t’aie  the  unsteady  pressure  on  an  osct.iatiss  re::aaiu,u: 

,..  ...  Ati-/subso  me  fiow,  :  b  e  seneralizcd-doublet-lattice- 
...  ,  .  j  - .  7  he  resuits  show  that  the  locally  linearized  panel 

..  ...yew  n-i curate  than  tlic  usual  linearized  methods.  But  further 

«. torts  are  sec-dcd  so  as  to  obtain  *  numerical  method  aTatiabie  Lr 
:  ii  luiatmq  tne  tran-.  supersonic  flows. 


PRESSURE  INTENSITY  DISTRIBUTION  COMPUTATION  FOR 
A  FLEXIBLE  WING  AT  TRANSONIC  SPEED 


Shen  Kewyang,  Zhuang  Xihua 

(Shanghai  Aviation  Industry  Corporation) 

ABSTRACT 

This  document  utilizes  an  iterative  problem  u 

solving  method  to  calculate  pressure  intensity  dis¬ 
tributions  of  a  flexible  wing  at  transonic  speeds. 

A  second  order  approximation  method  is  used  to  derive 
aerodynamic  force  for- small  perturbation  potential 
flow  at  transonic  speed;  a  one-dimensional  single 
beam  theory  or  a  two-dimensional  surface  structure 
matrix  is  used  to  calculate  a  flexion  coefficient  of 
effect.  All  results  shown  are  for  the  calculations 
of  a  M6  flexible  wing.  The  calculation  methods  for 
the  iteration  number  of  the  flexible  wing  and  the 
rigid  wing  are  the  same.  Flexion  change  for  wing 
surface  aerodynamic  force  load  distribution  and 
shock  wave  intensity  is  clearly  influenced. 

I.  Foreword 

As  all  of  us  know,  aircraft  wing  flexion  change  has  a  very 
important  influence  on  aerodynamic  force  load  distributions.  In 
transonic  speed  computations,  because  of  the  nonlinear  character¬ 
istics  of  the  fundamental  equation,  a  solution  for  aerodynamic 
force  is  obtained  through  an  iterative  method.  As  a  result  of 
this,  the  basic  line  of  thought  for  obtaining  a  solution  for 
transonic  speed  wing  flexion  aerodynamic  force  was  as  follows: 
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during  the  course  of  obtaining  an  iterative  solution  for  aero¬ 
dynamic  force,  an  interruption  in  the  sequenced  non-convergent 
aerodynamic  load  calculations  for  structural  change  occurred. 
Moreover,  the  quantitative  change  computed  during  the  course  of 
calculation  of  iteration  for  aerodynamic  force  showed  a  direct 
convergence.  According  to  this  reasoning,  Chipman,  et  al.,  [1] 
used  the  low  relaxation  method  of  iteration  to  create  an  equation 
for  computation  of  structural  change  at  transonic  speed  small 
perturbation  and  obtained  convergence  for  a  flexible  wing's  aero¬ 
dynamic  force  at  transonic  speed.  Whitlow,  et  al.,  [2]  took 
this  method  and  created  a  generalized  applied  equation  for  deri¬ 
vation  of  velocity  potential  at  transonic  speed. 

This  paper  utilized  the  transonic  speed  derivation  of  higher 
order  equation  (TSDH)  to  solve  for  aerodynamic  force.  Structural 
change  is  calculated  using  the  one-dimensional  simple  beam  model 
or  the  two-dimensional  surface  matrix  for  calculation  of  a  flexion 
coefficient  of  effect.  All  results  shown  herein  are  for  calcula¬ 
tion  of  an  M6  flexible  wing  and  they  have  gone  through  approxi¬ 
mately  100  iterative  derivation  trials  for  the  obtained  conver¬ 
gence  of  structural  change  and  pressure  intensity  distributions. 

At  the  same  time,  calculations  revealed  flexion  change  at  tran¬ 
sonic  speed  was  clearly  influenced  by  the  aerodynamic  forces. 


II.  Calculation  methodology 

After  designating  the  influx  Mach  number  as  Mm  and  designating 
velocity  pressure  as  q,  as  well  as  specification  of  the  wing 
surface  conditions,  the  aerodynamic  force  distribution  (l)  in 
small  perturbation  theory  is  merely  a  function  of  sectional  over¬ 
lap  angle  distribution  Cct)  and  sectional  slope  {dywfdx},.,  ,  so 

{!}=/({/>})  (1) 

In  the  equation  P  is  substituted  for  a  or  (dyjdx),.,  ,  and  here 
lower  point  u,l  is  separated  to  indicate  upper  and  lower  wing 
surfaces.  P  can  only  be  expressed  as 


(2) 


In  this  equation,  PR  expressed  the  values  of  rigid  characteris¬ 
tics  and  P£  expresses  the  contributions  of  structural  change. 

(PE>  is  merely  a  function  of  aerodynamic  force  (l) ,  so 

(3) 

An  iteration  method  can  be  used  to  solve  this  group  of  nonlinear 
equations.  This  paper  uses  the  low  relaxation  method  of  iteration 
for  calculation  of  structural  change,  so 

{PtY~“  -{PzY* +oLg<.{L)'") -{P;} '•']  (4) 

In  the  equation  (n)  and  (n+1)  distinguish  the  various  iteration 
values  for  n  trials  and  n+1  trials;  oi  is  a  factor  of  the  low 
relaxation  method  and  is  set  at  to  =  0.75. 


1.  Model  for  aerodynamic  force 


The  fundamental  equation  for  the  calculation  of  aerodynamic 
force  at  transonic  speed  is  the  TSDH  equation 

—  23/i(l 

The  corresponding  wing  surface  boundary  conditions  are  expressed 

**.<*.  ±0.r>-/,{(-fk-)#  |Cco*a-+n».(x,±0..-)3  -sinaj,  :***«<*>]  ^ 

In  the  equation,  e  is  a  parameter  for  small  perturbation,  e  =  «5*'7!/„, 
5  is  mean  correlational  thickness  of  the  wing.  1  is  the  small 
perturbation  velocity  potential,  y  is  the  comparative  specific 
heat,  f  is  an  adjustment  coefficient,  1  time) 

C 

(TN:  here  time  may  refer  to  measured  unit  or  trial) .  In  the 

obtuse  forward  peripheral  area  of  x  =  xLS(z)  ,  boundary  conditions  are 


(7) 


i 


In  this  equation,  -lu  is  the  reverse  angle  of  the  forward 
periphery.  The  obtuse  forward  peripheral  potential  velocity 
equation  is 
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i  -1-i  -*F 


—  2t§  au  ita.) 

\  a*  /  L  a*-  J 

/»•* 


ar.d  the  equation 


By  utilizing  all  of  the  difference  schemes  provided  in 
bibliographical  reference  [3]  for  the  equations  (5)-(8),  dis¬ 
cretization  for  the  difference  equation  group  may  occur. 
Furthermore,  after  using  the  scarce  density  mesh  substitution 
iteration  method  to  obtain  a  result  for  velocity  potential  $ 
convergence,  the  wing  surface  pressure  intensity  coefficient  C 

Sr 

can  be  used  for  calculation  in  the  formula  below: 

r 

+  +  -l]  (10) 

2 .  Structural  model 

(1)  One-dimensional  simple  beam  model  for  a  long,  tapered 
wing.  We  can  use  the  one-dimensional  simple  beam  (flexion  axis) 
theory  to  calculate  a  sectional  overlap  angle  for  elasticity  [4] 

f*  m<  *  ,  fj  )A/(  5  )</ 5  f*  /( z  ,  5  )T  i  5  )d  ij  mi 

a*(-’)  =  J.  -  Yl(f) +1  -5T<T> 

In  the  equation,  ar(  z )  is  the  sectional  overlap  angle  for  elas- 

A 

ticity  in  z.  *(:p5>  is  the  unit  moment  of  force  produced  for  the 
curve  distribution.  M'-.t})  is  the  aerodynamic  force  load  and 
inertial  load  distribution  rL)  along  the  curve.  £/(>}>  is  the  distri 
bution  of  cu^.e  rigidity  along  the  edge  of  the  flexion  axis. 

is  the  torsional  distribution  produced  by  the  unit  moment 

A 

of  force  within  z.  T(rj )  is  the  aerodynamic  load  (l)  and  the 
inertial  load's  torsion  distribution .  (?/( fj )  is  the  distribution 
of  torsional  rigidity  along  the  flexion  axis.  As  is  outlined 
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in  Figure  1(a),  lift  (l)  for  each  of  the  wing  sectional  plate 
components  perpendicular  to  the  z  axis  can  be  obtained  with 
integral  equation  (10).  Next,  by  using  equation  (11),  we  can 
calculate  each  plate's  flexion  overlap  angle .'{«,) , 


io.tr 


(a)  one-dimensional  simple 
beam  (flexion  axis)  model. 


(b)  two-dimensional  surface 
structural  model. 


Figure  1.  Structural  models  of  aircraft  wing 


(2)  Two-dimensional  surface  structural  model.  As  is  out¬ 
lined  in  Figure  1(b),  finite  trapezoidal  elements  are  divided 
into  specific  pressure  regions  on  the  wing  surface.  Through 
theoretical  calculations  or  experimentation,  a  coefficient  of 
effect  EC”]  for  flexion  is  obtained  for  each  surface  element. 
Hence,  the  quantitative  form  transformation  is  represented  as 

=  (12) 

In  the  equation,  (f’-  is  the  usable  amassed  load  of  the  form 

{F)  =  {'*CP-f:'>s}q  (13) 


In  this  equation,  is  obtained  from  equation  (10)  for  aero- 


fj  is  used 


dynamic  force  differences  (the  upper  orientation) 
for  the  form's  inertial  load  density  (the  lower  orientation), 
s  is  the  surface  of  the  finite  element.  Thus,  the  arising 
changes  of  the  wing  structure,  partitioned  into  upper  and  lower 
slope  surface  areas  are 

(14) 

(  Qy.  \  \  _  ay- 

\  OX  ■  £■)  .  ;  Ox 


{(- 


In  this  manner,  through  use  of  equation  (11)  or  equation 


(14),  we  can  obtain  a  value  for  rPT?}  for  use  in  equation  (2). 


In  computation,  for  each  iteration  of  aerodynamic  force  m  times 
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(for  loose  mesh,  m  =  10;  for  packed  mesh,  m  =  5) ,  it  is  con¬ 
venient  to  use  equation  (11)  or  equation  (14)  to  figure  f P  } 
initially.  After  obtaining  a  value  for  (p)  in  equation  (4) 
using  the  low  relaxation  method,  the  value  can  then  be  substi¬ 
tuted  in  equation  (6)  to  calculate  the  next  m  iterations  of 
aerodynamic  force  until  the  aerodynamic  force  and  the  value  (Pg) 
both  converge. 

Ill .  Calculation  results  and  discussion 

In  order  to  test  the  feasibility  of  this  method  ,  we  selected 
known  finite  rigidity  results  for  an  M6  wing  (M^  =  0.84,  a  = 

3.06°)  (3]  as  a  calculation  model.  We  used  a  1/2  wing  spread 

2 

value  of  7.41  meters,  influx  velocity  pressure  q  was  5128  kg/m  . 
The  flexion  axis  was  at  a  50%  tension  field,  the  curve  rigidity 
distribution  El  is  shown  in  Figure  2.  Torsional  rigidity  was 
not  figured  in  the  calculation.  In  regards  to  the  two-dimensional 
structural  model,  the  flexion  coefficient  of  effect  is  calculated 
using  the  equation  below; 

C"< 

In  this  integral,  che  -pper  limit  z  is  fj  > :  when  z  =  z,  and  when 

/ 

l  <  z  it  is  z  = 

The  following  are  calculation  conditions  for  the  model: 

(1)  Flexible  wing  M6E1  (a  =  3.06°) — =0.84,  a  =  3.06°. 

In  regards  to  the  influence  of  pressure  intensity  distributions, 
we  must  understand  the  overlap  angle  does  not  change  as  flexion 
change  takes  place.  El  expresses  use  of  the  simple  beam  theory 
in  equation  (11)  for  the  calculation  of  flexion  overlap  angle. 

(2)  Flexible  wing  M6E1  (a  =  5.2°)  — =0.34,  a  =  5.2°. 

This  is  the  flexible  wing's  lift  and  the  origin  for  rigid  wing 

M6  comparisons.  In  regards  to  the  influence  of  pressure  intensity 
distributions,  we  must  understand  that  lift  does  not  change  as 
flexion  change  takes  place. 


Ctj  ~  1  >  <  *  -Jj. 


TU) 


(3)  Flexible  wing  M6E2  (a  =  3.06°) — Mm  =  0.84,  a  *  3.06°. 

E2  expresses  a  two-dimensional  flexion  matrix  used  in  equation  (2) 
to  calculate  structural  change.  This  calculation  process  uses 
the  data  of  the  two-dimensional  structural  model  to  postulate 
reasonable  possibilities  for  our  method.  Furthermore,  the  cal¬ 
culated  results  are  comparable  to  those  for  M6E1. 

The  results  of  these  calculations  are  illustrated  in  Figures 
3-6.  The  iterative  calculation  process  converges  on  the  objective 
£  =  MAX -33/;;*)  |<  10'3  •  Bibliographical  reference  [5]  further 
expands  upon  the  results  of  this  experimental  objective,  and  it 
provides  wind  tunnel  data  to  coincide  with  it. 

Figure  3(a)  expresses  the  wing  flexion  overlap  angle  a£  for 
M6E1  under  conditions  of  change  as  an  iteration  number  over  n 
trials.  As  you  can  see,  stability  increases  with  an  increase  in 
the  change  in  the  iteration  number.  The  change  for  M6E2  is  also 
likewise.  Figure  3(b)  expresses  the  wing  flexion  overlap  angle 
a£  for  changes  that  occur  in  n .  As  the  figure  clearly  shows, 
a  two-dimensional  structural  model's  flexion  axis  deflection  curve 
is  obtained  by  the  flexion  overlap  angle  and  is  close  to  the 
results  for  M6E1,  so  we  have 

a*  -  —  s  i  n  A 

a  z 

in  the  formulae,  A  is  the  flexion  axis  reverse  sweep  angle,  dy-.dz 
is  the  slope  change  for  the  flexion  axis  along  the  flexion  axis 
deflection  surface. 

Figures  4  and  5  differentiate  results  of  chord  direction  for 
pressure  intensity  distributions,  sectional  lift  coefficient  C^, 
and  the  sectional  pressure  point  x^  for  three  aforementioned 
models.  Furthermore,  these  results  are  compared  to  the  results 
for  M6 .  As  you  can  see  in  the  figure,  aerodynamic  force  flexion, 
aside  from  the  serious  influences  of  lift  distribution,  is  clearly 
influenced  by  pressure  intensity  distributions  and  shock  wave 
intensity.  At  the  same  time,  the  closeness  of  results  for  M6E1 
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(ot  =  3.06°'  and  M6E2  (tc  =  3.06°)  remained  true  for  all  mathematical 
theory  me  the  is  tried  in  experiments  for  a  ev;c  -dimensional  struc¬ 
tural  mode  _  J-* 


a:  NsEi 

Figure  2.  M6S2  curve  rigid¬ 
ity  coin.'iiir.g  to  matrix 
changes 

1 — e1*10  3  kg/m;  2 — meters 
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Figure  3.  Structural  change 
coinciding  with  wing  change 
and  iteration  number. 

1 — loose  mesh;  2 — packed  mesh; 
3--increased  density;  4 — even 
more  increased  density; 

5--M6E2 (a=3.06°)  via  flexion 
axis  deflection  curve  of  ob¬ 
tained  flexion  overlap  angle. 


F  gure  4,  Comparison  of 
pressure  intensity  computa¬ 
tion  results. 
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Figure  5.  Comparisons  of 
secoional  lift  coefficients 
and  relative  pressure. 


TA3LE  1.  Aerodynamic  force 
derivative  comparisons 


1 

z  — ,4 

,  Ms  V6E1/M6E2  K  =  ~-=-3_ 

J  c;,:7* 

0.C431  3.  *70 

^  m’r.l/Z. 

-0.333-J  C.537 

1 — aerodynamic  derivative;  z-- 
k=flexion  value/rigiditv  value; 
3 — 1/ degrees 

TABLE  2.  Iterative  trial 
comparisons  (e=10“M 
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l — conditions 


Finally,  Table  1  shows  how  aerodynamic  force  flexion  is 

A  A 

influenced  by  linear  slope  of  lift  and  derivative  Mz  for  the 
pitching  moment.  The  table  shows  a  significant  difference  for 
M6E1  and  M6E2.  Table  2  shows  comparisons  for  the  iterative 
number . 


As  you  can  see,  the  iteration  numbers  for  a  flexible  wing 
and  a  rigid  wing  are  nearly  equal.  Resultingly,  the  flexible 
wing  transonic  speed  computations  take  about  the  same  amount  of 
time  to  figure  as  for  a  rigid  wing,  using  less  than  one  hour  of 
CPU  time  on  an  IBM  4341  computer. 


IV.  Conclusions 

For  the  nonlinear  nature  of  near  transonic  speed  aerodynamic 
force  computations,  this  paper  in  the  course  of  iterative  deriva¬ 
tion  for  aerodynamic  force,  interrupts  the  computation  of  the 
aerodynamic  force  based  on  its  previously  obtained  values  for 
structural  change,  then  these  values  for  quantitative  change  are 
placed  into  the  iteration  process  for  aerodynamic  fcrce  until  they 


converge.  For  the  scope  of  M6  flexible  wing  calculations  shown 
herein,  this  computation  method  proved  quite  effective  in  find¬ 
ing  solutions  for  a  flexible  wing's  quantitative  change  and  aero¬ 
dynamic  force  distribution  convergence.  Because  of  this,  the 
method  cited  herein  can  be  used  to  predict  the  aerodynamic  force 
and  flexion  for  the  lift  plane  at  transonic  speed  as  well  as  to 
predict  aerodynamic  flexion  cuts  in  the  lift  surface  design. 

This  research  work  was  supported  by  the  Engineering  Division 
of  the  Ministry  of  Aviation  Industry  and  we  express  our  gratitude 
for  their  assistance. 
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TRANSONIC  PRESSURE  DISTRIBUTION  COMPUTATIONS 
OF  A  FLEXIBLE  WING  - 

Shen  Keyang,  £hang  X  i  h  u  a 
■  Shanghai  Aviation  Industry  Corjcrsf.cn. 

Abstract 

Aa  iteration  method  is  used  tor  computing  the  transonic  pressure 
distribution  on  a  flexible  wing.  The  aerodynamic  force  is  solved 
with  the  second  order  approximation  method  for  transonic  small 
disturbance  potential  flow.  The  structural  deflection  is  computed  using 
the  one-dimensional  simple  beam  theory  or  a  two-dimensional  flexibility 
matrix  method.  The  typical  computations  for  Ms  flexible  wing  indicate 
that  the  iteration  number  for  flexible  wing  computation  is  only  slightly 
more  than  that  for  rigid  one,  and  that  the  influences  of  the  structural 
deflection  on  the  airload  and  shock  strength  are  important. 
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ON  RELAXATION  OF  TRANSONIC  VELOCITY  ZERO  LIFT  AIRFOIL 
POTENTIAL  FLOW  AND  CONVERGENCE  OF  TWO  SELF-CORRECTING 
WIND  TUNNELS 


Liu  Xueding,  Luo  Shijun 
(Northwestern  Polytechnical  University) 

ABSTRACT 

This  document  utilizes  a  transverse  small  per¬ 
turbation  velocity  potential  equation  or  a  longitu¬ 
dinal  large  perturbation  velocity  potential  equation 
to  calculate  the  potential  flow  for  a  zero  lift 
airfoil  at  transonic  velocity.  In  numeric  line 
relaxation  testing  [1],  when  the  difference  equa¬ 
tion  uses  a  simple  iteration  methodology  and  when  the 
4>xx  difference-  equation  uses  an  improved  iteration 
methodology,  the  characteristics  of  stability  compare 
favorably.  This  conclusion  coincides  with  the 
theoretical  analysis  of  linearization  cited  in  biblio¬ 
graphical  reference  [2]. 

This  document  uses  a  mixed  difference  method  of 
numeric  simulation  to  verify  the  static  pressure 
along  the  wing  surface  in  regards  to  a  two  control 
surface  for  convergence  of  zero  lift  airfoils  in  a 
self-correcting  wind  tunnel  at  transonic  velocity. 

For  the  former  scheme,  a  NACA  0012  airfoil  is  used, 

M  =  0.9  and  a  RAE  104  airfoils  is  used,  =  0.8. 

In  regards  to  the  former  scheme,  a  NACA  0012  airfoil 
is  used,  Mw  =  0.72,  0.8.  When  the  stream  angle  is  at 
zero  and  wind  tunnel  height  is  comparable  to  wing 
chord,  we  can  nearly  get  a  convergence  of  free  flow 
with  no  wind  tunnel  wall  interference  after  about 
three  hours. 


This  document  received  July  19,  1983.  Revisions  received 
January  5,  1984. 
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When  using  the  mixed  difference  line  relaxation  method  to 
solve  transonic  potential  flow  equations,  we  ran  into  a  problem 
with  convergence.  The  majority  of  computation  experiments  [1] 
revealed  the  following:  during  calculations  involving  super¬ 
critical  flow,  the  iteration  method  caused  a  problem  that 
Influenced  the  stability  for  convergence.  This  document  discusses 
the  problem  and  the  experimental  test  results  which  turned  out  to 
be  significant  for  it. 

Bibliographical  reference  [3]  cites  a  self-correcting  method 
that  uses  two  airflow  parameters  and  <f>y  for  a  one  control 
surface.  Numeric  simulation  proved  that  when  supercritical  flow 
manifests  itself  within  the  field  of  flow  for  the  NACA  65A003 
airfoils,  the  self-correcting  wind  tunnel  converged.  Considering 
the  difficulties  of  obtaining  a  measurement  of  <J>  in  wind  tunnel 
experiments,  this  document  proposes  that  only  two  measurements, 
for  static  pressure  and  <t>x,  need  to  be  done  in  the  self-correcting 
wind  tunnel  scheme.  The  convergence  was  proven  in  this  manner 
by  utilising  numeric  simulation. 


I .  An  iteration  method  for  the  influences  of  convergence 


Consider  a  normal  surface  potential  flow  at  transonic  speed. 
The  hypothetical  orientation  perturbation  velocity  y,  as  compared 
to  the  last  crossed  oerturbation  velocity  of  q  is  much  less. 
However,  the  x  orientation  perturbation  velocity  has  the  possi¬ 
bility  of  not  being  less.  Perturbation  velocity  can  be  adequately 
represented  in  a  partial  differential  equation  as 

<1  — (1.1) 

and  in  the  formula 


1-A/1 


(1.2) 


y  is  the  airflow's  adiabatic  index  number,  Mm  is  the  last  crossed 
Mach  number  for  perturbation  airflow. 


Equation  (1.1)  is  a  mixed  form  nonlinear  equation.  Due  to 

the  multi-faceted  nonlinear  characteristics  of  the  group  of 

equations  as  seen  in  (1.1),  we  can  use  a  line  relaxation  iteration  /35 

method  to  linearize  a  solution.  As  the  edge  of  the  y  axis  orients 

towards  line  relaxation,  an  iteration  equation  is  required  for 

2 

the  velocity  discriminant  1-M  ,  $  coefficient  and  the  various 

aA 

$  values  of  velocity  potential  apart  from  the  relaxation  line. 

AX 

According  to  page  62  of  f4],  an  accurate  difference  scheme 
for  velocity  discriminants  can,  on  the  one  hand,  use  a  simple 
iteration,  but  for  supercritical  airflow,  this  can  lead  directly 
to  a  divergence  of  vibration.  This  can  be  organized  into  several 
iterative  methods  as  cited  below: 

simple  iteration,  $  simple  iteration  (P  =  0) 

X  XX 

,  .  .  <p  improved  iteration,  <J>  improved  iteration  (P  *  1) 

velocity  discriminant  x  xx 

1-M*simple  iteration  d>x  simple  iteration,  $xx  improved  iteration  (P  *  2) 

<f>  improved  iteration,  $  simple  iteration  (P  *  3) 

X  XX 

Numeric  experimentation  [1]  revealed  the  following  informa¬ 
tion:  for  M  <0.775  free  flow,  for  use  of  the  P  =  0,  1  and  2 

iterative  methods,  the  calculations  approximated  convergence.  The 
only  difference  was  in  the  speed  at  which  convergence  occurred. 

For  Mw  ^_0,8  free  flow,  the  supposed  local  supersonic  region  was 
relatively  large,  and  an  intense  shock  wave  manifested.  Using  the 
two  iterative  methods  P  =  0,  1,  it  was  extremely  difficult  for  a 
convergence  to  occur  in  computation.  However,  using  the  iteration 
method  with  P  =  2,  a  calculation  result  that  systematically  led 
to  success  with  was  obtained.  The  greatest  value  for  was 
about  0.925.  With  a  39x23  mesh,  the  Ximen  7760  computer  took 
about  10  minutes  to  perform  the  calculations.  This  goes  to  say 
that  the  iteration  method  for  P  =  2  are  one  way  to  use  computations 
of  supercritical  flow  to  make  comparisons  of  stability.  The 
theory  and  results  of  this  discussion  is  analyzed  in  12].  Based 
on  [2],  the  simple  iteration  method  for  <p  localized  supersonic 
points  is  unstable.  For  the  condition  of  P  =  3,  this  paper  per¬ 
formed  no  experiments  or  calculations. 
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II .  Self-correcting  wind  tunnel  convergence  for  zero  lift 
symmetric  airfoils  in  view  of  static  pressure  at  transonic 
speed  for  two  control  surfaces 


This  document  uses  the  same  airflow  parameters  <J>  static 

X 

pressure  for  two  control  surfaces  in  the  self-correcting  wind 
tunnel.  The  proofs  given  are  in  regards  to  the  NACA  0012  airfoil  RAE  104 
When  supercritical  flow  manifests  within  the  field  of  flow,  there 
is  convergence  within  the  self-correcting  wind  funnel.  A  mathe¬ 
matical  model  is  depicted  in  Figure  1  below: 


Key :  a — external  field  of  flow  boundary;  E — foward  field  of  flow 
boundary;  C — aft  field  of  flow  boundary;  D — region  IT;  E — wind 
tunnel;  f — second  control  surface;  0 — common  region;  H — first 
control  surface;  I — region  1;  j — main;  X — wind  tunnel  axis; 

L — stream  angle;  M — main 


By  taking  the  first  and  second  control  surfaces  as  boundaries 
the  fields  of  flow  are  broken  into  regions  I  and  II  and  a  common 
region  III  or  three  components  altogether.  Region  I  plus  the 
common  region  makes  up  the  interior  field  of  flow.  Region  II 
plus  the  common  region  makes  up  the  exterior  field  of  flow. 

In  theory,  region  II  should  approach  a  boundary less  state.  The 
first  and  second  control  surface's  parameters  of  flow  are 
separated  into  and  •  Lower  coordinate  regions  E  and  T 
are  separated  to  express  an  external  field  of  flow  and  an  internal 
field  of  flow.  N  is  the  iteration  number  for  the  self-correcting 
wind  tunnel. 


We  start  our  calculations  from  the  solution  for  the  internal 


field  of  flow.  The  external  boundary  conditions  for  the  initial 


(0) 


field  of  flow  are  =  0  w^ich  is  f°r  tbe  boundary  conditions 

at  the  mouth  of  the  wind  tunnel.  In  computations  for  the  inter¬ 


nal  field  of  flow  for  the  first  control  surface,  the  distribution 


(N) 


of  perturbation  velocity  for  the  x  orientation  is  $xlT  and  this 
creates  the  internal  boundary  conditions  for  the  external  field 
of  flow,  so 


1 1  ^PuT 


In  computations  for  the  external  field  of  flow  for  the  second  con- 

ribution  perti 
As  an  example 

So  we  have  -  a ~K)i>i Vr 


trol  surface,  the  distribution  perturbation  velocity  for  the  x 
*  (N) 
vx2E* 


orientation  is 


Consecutive  calculations  are  made  for  the  internal  field  of  flow's 
external  boundary  conditions.  Within  the  self-correcting  wind 
tunnel,  relaxation  factor  was  preset  to  the  value  of  K  =  0.5  for 
purposes  of  this  paper.  For  iteration,  this  factor  led  to  a  direct 
convergence . 


The  flow  of  calculation  is  depicted  in  Figure  2. 


Figure  2.  Flow  of  numeric 
simulation  calculations  for 
two  surfaces  self-correcting 
wind  tunnel  static  pressure. 
Key:  A — selected;  2 — (initial 

conditions)  ;  C — internal  field  of 
flow  difference  iteration  calcula¬ 
tions;  D — external  field  of  flow 
difference  iteration  calculations; 
E — calculation  of  airfoil  surface 
pressure  coefficient;  F — true; 
-false 


The  boundary  conditions  for  the  internal  field  of  flow  are 
expressed  below: 


forward  field  of  flow  boundary  =  3 

aft  field  of  flow  boundary  <p  =  3 

(0)  x 

second  control  surface  =  0  (initial  conditions) 


Wind  tunnel  axis  line 


$  mi?'  i  - 

<  1  -  K  )<r  ,'ir 


(2.1) 


air  foil  exterior  <J>  =0 

y 

point  on  forward  edge 

top  of  airfoil  when  stream  angle  is  an  zero 

ax 

The  boundary  conditions  for  the  external  field  of  flow  are 
cited  below: 

forward  field  of  flow  boundary  <J>  =  C 
aft  field  of  flow  boundary  6^  =  0 

upper  field  of  flow  boundary  $  =  0 

first  control  surface  =  ^xlT  (&■"& 


The  aforementioned  boundary  conditions  for  (2.1)  and  (2.2) 
may  be  altered  to  take  the  form  below: 

The  internal  field  of  flow's  x  orientation  mesh  interval  is  Ax, 


<P!r,,-Krir^u-K)<F!r  (2#la) 

and  *'*  V,T  (2.2a) 

It  is  much  more  convenient  to  use  (2.1a)  and  (2.2a)  instead  of 
(2.1)  and  (2.2.) . 


ff,  <  V>  Cjp 
r  I  t  ““  *F  r  r 


Calculation  results  are  depicted  in  Figures  3-6.  For  the 
NACA  0012  airfoil,  0.9  and  P=2  iteration  method  was  used. 

For  the  RAE  104  airfoil,  Mm  =  0.8,  and  the  P=1  iteration  method 
was  used.  The  internal  field  of  flow  and  the  external  field  of 
flow  were  separated  for  relaxation  iteration.  After  two  itera¬ 
tion  trials,  convergence  on  the  objective  occurred.  The  minimum 
absolute  values  for  perturbation  velocity  potential  of  each  node 
was  10  ^ ,  and  the  maximum  values  were  less  than  10  \  At  this 
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Figure  3.  Airfoil  pressure 
distribution  comparisons  for 
the  NACA  0012  in  free  flow  and 
self-correcting  wind  tunnel. 

A — free  flow;  B--self-correction 
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Figure  4.  Comparisons  of  control  surface  pressure 
distribution  for  NACA  0012  airfoil,  *  0.9,  u  =  C. 3 
in  free  flow  and  self-correcting  wind  tunnel. 

A — free  flow;  B — self-correcting  wind  tunnel  far  wall  control 
surface;  C — self-correcting  wind  tunnel  near  wall  control 
surface 


point,  after  the  separation  convergence  for  the  external  field  of 
flow,  we  made  a  comparison  of  it  with  the  nodal  points  of  the 

second  control  surface  for  perturbation  velocity  within  the  flow. 

-4 

When  the  greatest  absolute  difference  value  obtained  was  10  ,  it 

was  recognized  as  a  convergence  to  no  wind  wall  disturbance  within 
the  seif -correcting  wind  tunnel. 

Calculation,  results  for  free  flow  and  self-correcting  wind 
tunnel  static  pressure  distributions  for  both  types  of  airfoils 
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Figure  5.  RAE  104  airfoil  pressure 
distributions  in  free  flow  and 
self-correcting  wind  tunnel  (M  = 

0 . 8 ,  a  =  0  ° )  . 

A — free  flow;  B--self-correction 


are  depicted  in  Figures  3  and  5.  The  results  obtained  for  both 
calculation  methods  are  similar. 

Calculation  results  for  free  flow  and  self-correcting  wind 
tunnel  pressure  distributions  for  two  control  surfaces  are  com¬ 
pared  in  Figures  4  and  6.  The  results  obtained  for  both  of 
these  methods  is  nearly  identical. 


c,  k« 


Figurs  6.  Comparisons  of  control  surzace  pressure 
distributions  for  RAE  104,  M  =0.8,  a  =  0°,  in  free 
flow  end  self-correcting  wind  tunnel. 

A — free  flow;  3--self-correcting  wind  tunnel  far  wall  control 
surface;  C — self-correcting  wind  tunnel  near  wall  control 
surface 


In  consideration  of  measurement  and  discernment  of  two 
control  surface  static  pressure  distributions,  there  are  certain 
difficulties.  This  article  proposes  a  different  scheme  for  self 
correcting  wind  tunnels.  A  mathematical  model  is  depicted  in 
Figure  7. 
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Figure  7.  Principles  of  self-ccrr  _  -  .ir.g  wind  tunnel 
as  applied  ~o  control  surface  ar.i  i-.-fcil  static 
pressure . 

A— upper  boundary  of  external  field  of  flow;  B — region  II; 

C — tunnel  wall;  D — control  surface;  E — region  I;  F — wind  tunnel 
axis;  G — stream  angle 


c.r 


Figure  8.  Comparison  of  airfoil  pressure;  distributions  in 
free  flow  anc  self-correcting  wind  tunnel  for  NACA  0012, 
=0.72,0=0°. 

A — free  flow;  B — self-correction 


By  taking  the  control  surface  as  a  boundary,  we  can  separate  the 
field  of  flow  into  two  components,  region  I  and  region  II. 

Region  I  is  the  internal  field  of  flow.  Region  I  plus  region  II 
comprise  the  external  field  of  flow.  In  theory,  region  II 
nears  an  unlimited  or  unbounded  state.  The  atmospheric  flow 
parameters  for  the  control  surface  and  the  airfoil  are  divided 
into  V,K  and  .  Lower  boundaries  E  and  T  express  external 

flow  and  internal  flow  fields.  N  is  the  self-correcting  wind 
tunnel  iteration  number. 
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Figure  9.  Comparison  of  air 
foil  pressure  distributions 
in  free  flow  and  self-correct¬ 
ing  wind  tunnel  for  NACA  0012, 
=  0.72,  a  =  0°. 

A — free  flow; 

B — self-correcting  wind  tunnel 


*  ^  .*  x/e  i 

Figure  10.  Comparison  of 
airfoil  pressure  distributions 
in  free  flow  and  self-correct¬ 
ing  wind  tunnel  for  NACA  0012, 
Mm  =  0.8,  a  =  0°. 

A — free  flow 
B — self -correction 


Starting  with  calculation  of  a  solution  for  the  internal 
field  of  flow,  the  initial  internal  field  of  flow’s  external 
boundary  conditions  are  *  0.  From  the  internal  field  of 

flow,  we  calculate  the  x  orientation  perturbation  velocity  dis¬ 
tributions  for  the  airfoil  as  By  using  as  the  internal 

boundary  conditions  for  the  external  field  of  flow,  we  have 


From  within  the  external  field  of  flow,  we  can  calculate  the  con¬ 
trol  surface's  x  orientation  perturbation  velocity  distribution 


<t> 


(N)  . 


xKE  * 

For  example 

which  yields  ?*!?"  mKvltt+  (1— 

which  becomes  the  internal  field  of  flow's  external  boundary 
conditions  on  successive  calculations.  Herein,  K  is  the  relaxa¬ 
tion  factor  for  the  self-correcting  wind  tunnel.  This  article 
set  the  value  of  K  =  0.5.  For  iteration,  it  led  to  a  direct 
convergence.  The  calculation  model  is  depicted  in  Figure  2. 


The  internal  field  of  flow  boundary  conditions  is  expressed 

below; 
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self-correcting  wind  tunnel  is  in  Section  2.  For  the  NACA  0012 
airfoil,  =  0.72,  0.8,  a  =  0°.  Both  calculations  for  free  flow 
and  self-correcting  wind  tunnel  are  compared  in  Figures  8  and  10. 

For  both  computations,  the  results  are  similar-  Both  cal¬ 
culations  for  control  surface  pressure  distributions  in  free 
flow  and  self-correcting  wind  tunnel  are  compared  in  Figures  9 
and  11.  The  results  for  both  calculation  methods  are  nearly 
identical . 

IV .  Conclusion 

Based  on  numeric  experimentation,  the  findings  of  this 
document  show  improved  stability  characteristics  for  the  longi¬ 
tudinal  large  perturbation  and  transverse  small  perturbation 
transonic  flow  around  a  symmetric  airfoil  in  a  zero  stream  angle 
and  in  linear  relaxation  where  the  difference  formula  for  <!> 
uses  improved  iteration  and  the  difference  formulas  for  4>x  uses 
a  simple  iteration. 

In  regards  to  transonic  velocity  zero  lift  symmetric  air¬ 
foils  and  self-correcting  wind  tunnels,  this  document  proposes 
using  the  self-correcting  method  for  static  pressure  for  a  two- 
control  surface  airfoil  and  the  self-correcting  method  for  static 
pressure  on  a  one-control  surface  airfoil.  In  this  way,  numeric 
simulation  can  prove  a  convergence. 

The  research  of  this  document  c~n  be  generalized  to  tri- 
elemental  flow  in  a  non-zero  lift  state. 
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ON  RELAXATION  OF  TRANSONIC  FLOWS  AROUND 
ZERO-LIFT  AIRFOIL  AND  CONVERGENCE  OF 
SELF-CORRECTING  WIND  TUNNELS 
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Liu  Xueding,  Luo  Shijun 
(Northwestern  Polytechnical  University)  •  - 

Abstract 

With  the  assumption  of  small  transverse  velocity  components,  the 
steady  traasonic  potential  flow  around  symmetric  airfoil  at  zero  angle  of 
attack  is  computed  by  the  mixed  difference  method.  After  some  numer¬ 
ical  experiments  on  the  stability  of  various  possible  schemes  of  iteration 
in  the  relaxation,  a  stable  scheme  is  found  and  used  to  verify  the 
convergence  of  two  kinds  of  transonic  seif-correcting  wind  tunnels 
which  are  based  on  the  pressure  distributions  along  (1)  two  control 
surfaces  and  (2)  one  control  surface  and  the  airfoil. 


OBSERVATION  OF  FLOW  DIAGRAMS  FOR  EXPANSION  SHOCKWAVE  TUBE 
Yu  Hongru,  Lin  Jianmin,  Li  Zhongfa  and  Gu  Jiahua 
1 .  FOREWORD 

The  calculation  methodology  for  uniform  shockwave  tubes  has 
already  developed  to  near  perfection.  In  order  to  obtain  an  intense 
shockwave,  we  often  make  use  of  a  "contracting"  shockwave  tube. 

Its  booster  section  is  about  the  same  as  a  booster  stage.  A  diaph¬ 
ragm  has  been  added  to  allow  the  nozzle  to  contract.  A  large 
booster  stage  is  only  required  in  research  on  low  density  gases 
or  in  low  intensity  shockwave  research.  In  order  to  reduce  the 
quantity  of  gas  used  in  the  booster,  the  "expansion"  shockwave 
tube  is  generally  selected  for  utilization.  Its  booster  section 
is  smaller  than  a  booster  stage,  and  expansion  through  the  nozzle 
is  similar  for  both  types.  In  regard  to  the  shockwave  tube  move¬ 
ment  at  the  slotted  boundary  region,  a  contraction-expansion 
nozzle  is  often  placed  in  the  lower  portion  of  the  diaphragm  to 
prevent  the  mutual  interference  from  the  rarefaction  wave  and  the 
reflex  shockwave  from  occurring  or  it  is  done  in  order  to  obtain 
a  polar  extreme  shockwave.  As  for  this  type  of  variational  sec¬ 
tional  shockwave  tube  flow,  Resler,  et  al.  II]  states  that  for  the 
contracting  shockwave  tube,  it  becomes  a  matter  of  sonic  or  super¬ 
sonic  flow  as  it  passes  through  the  lower  nozzle  area,  and  he  pro¬ 
vided  a  calculation  msthodology  for  shockwave  intensity.  Alpher 
and  White  [2]  took  the  methodology  provided  in  bibliographical 
reference  [1]  and  added  to  it,  making  several  improvements.  They 
advocate  use  of  a  diaphragm  with  a  nozzle  that  can  be  adjusted  for 
various  types  of  flow,  and  not  only  can  the  nozzle  be  set  for 
sonic  or  supersonic  flow,  but  furthermore,  the  nozzle  can  be  set 
for  subsonic  flow  as  well. 

If  the  transition  nozzle  is  of  the  expansion  type,  Roscizewcki 
13]  points  out:  when  the  nozzle  is  at  supersonic  flow,  it  exists 
in  two  different  and  distinct  types  of  flow  wave  diagrams — as  pre¬ 
expansion  in  the  lower  reaches  manifests  itself  as  a  rarefaction 
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wave,  and  as  post  expansion  manifests  itself  as  a  secondary  shock- 
wave.  The  two  aforementioned  theories  are  both  based  on  regulated 
flow  of  the  booster  gases.  This  is  because  it  is  not  suitable  to 
use  with  the  secondary  shockwave  flow  that  ensues.  Laderman  [4] 
provides  a  method  of  flow  calculation  for  the  secondary  shockwave 
that  comes  about. 


Russel  [5]  makes  a  diagnosis  for  expansion  shockwave  tube  flow 
The  theory  for  conditions  of  pre-expansion  predicts  that  the  rare¬ 
faction  wave  in  the  lower  reaches  of  the  nozzle  will  be  very  easily 
observed.  However,  the  theory  for  conditions  of  post  expansion 
secondary  shockwaves  in  the  lower  reaches  of  the  nozzle  were  not 
observed  as  predicted. 


The  contents  of  this  document  show  successive  experimental 
observation  results  for  expansion  shockwave  tube  flow  as  well  as 
the  corresponding  flow  diagrams  for  actual  conditions  of  flow  as 
they  exist. 


Figure  1.  Four  types  of  wave 
flow  systems  for  shockwave 
tubes 

a — pre-expansion  flow,  secondary 
shockwave  as  rarefaction  wave; 
b — post-expansion  flow,  secondary 
wave  as  shock  or  compression  wave; 
c — post-expansion  flow,  secondary 
wave  in  nozzle  interior;  d — sub¬ 
sonic  velocity  flow;  e — primary 
wave;  f — boundary  region;  g — sec¬ 
ondary  wave;  h — rarefaction  wave 
at  booster  interior 


2.  DIAGNOSTIC  METHOD 


The  expansion  shockwave  tube  exists  in  four  different  forms 
of  flow  conditions  (Figure  1)  . 


a.  When  the  primary  shockwave  and  the  Mach  number  exceed  a 
specified  critical  value,  a  rarefaction  wave  manifests  in  the  lower 
reaches  of  the  nozzle.  This  particular  type  of  flow  is  referred 


•'Utica.  c»aA««at 


ci^ua 


to  as  pre-expansion  flow.  If  a  pressure  blockage  occurs  within 
the  lower  reaches  of  the  nozzle  tube,  the  pressure  curvilinear 
conditions  and  formulas  are  as  expressed  in  Figure  2.  Following 
a  decrease  in  the  intensity  of  the  primary  shockwave,  the  rare¬ 
faction  wave  begins  its  decay  process.  At  the  time  when  the 
rarefaction  wave  has  equalized  to  that  of  the  Mach  wave,  we  refer 
to  this  moment  as  the  "greatest  disturbance"  of  expansion. 


Figure  5.  Realm  of  four 
types  of  expansion  shock 
wave  tube  flow  fields. 


1 —  pre-expansion  region; 

2 —  H  N  propellant; 

3 —  N  N  propellant; 

4 —  post-expansion  region; 

5 —  pre-expansion  region; 

6 —  post-expansion  region; 

7—  post-expansion  region; 

8 —  subsonic  area; 

9 —  subsonic,  region 


b.  When  the  primary  shockwave 
and  the  Mach  wave  decrease  to  "greatest 
disturbance",  the  primary  shockwave 
transitions  to  a  post  expansion  flow. 

If  Roscizewcki  foresaw  that  particular 
type  of  manifestation  occurring  as  a 
secondary  shockwave,  the  pressure  curve 
conditions  and  corresponding  wave  flow 
diagrams  would  be  as  depicted  in  Figure 
3a.  Due  to  the  disturbance  that  exists 
within  the  boundary  layer  of  actual 
flow,  it  is  clear  that  it  can  only 
exist  as  a  compression  shockwave.  At 
this  time,  the  curve  measurements  for 
pressure  conditions  and  the  correspond¬ 
ing  wave  flow  diagrams  would  be  as  are 
shown  in  Figure  3b. 


c.  As  the  intensity  of  the  primary  shockwave  successively 
decreases,  the  secondary  shockwave  or  compression  wave  moves  into 
the  interior  of  the  nozzle. 


d.  With  a  further  decrease  in  the  intensity  of  the  primary 
shockwave,  the  booster  gases  change  completely  and  become  subsonic 
flow. 


The  contraction  shockwave  can,  therefore,  only  exist  under 


flow  conditions  a  and  b 
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This  experiment  is  limited  to  observation  of  the  two  types  of 
flow  depicted  in  wave  diagrams  a  and  b.  From  a  comparison  of 
Figures  2  and  3,  we  can  see  the  following:  A  differentiation  of 
pre-expansion  flow  and  post-expansion  flow  based  on  whether  or  not 
the  aft  platform  is  positioned  higher  than  the  forward  platform 
(mounting  platform) .  The  result  of  a  post-expansion  flow  is  a 
secondary  wave  or  a  compression  wave  which  is  reliant  upon  the 
existing  pressure  curve  conditions.  If  we  can  overlook  the  influences 
of  the  boundary  layer  and  the  sensing  devices,  the  secondary  shock 
wave  has  a  steep  trailing  edge;  furthermore,  the  trailing  edge  of 
the  compression  wave  that  is  produced  has  a  definite,  specific 
width.  In  actual  measurements,  there  was  a  specific  dimension 
given  for  the  sensing  devices,  and  the  secondary  shockwave  passed 
through  the  sensing  device  using  up  a  specified  amount  of  time. 

This,  coupled  with  the  primary  shockwave  that  previously  occurred, 
creates  our  boundary  layer.  As  this  fork  or  division  occurs,  a 
certain  amount  of  time  is  required  for  it  to  pass  through  the 
sensing  devices.  As  a  result,  there  is  a  gap  of  a  specific  width 
in  the  measurements  taken  for  the  pressure  curve.  However,  for  a 
comprehensive  analysis,  the  two  widths  cannot  be  the  same.  We  must 
base  our  diagnosis  on  the  dimensions  of  one  of  these  widths. 

3.  EXPERIMENTAL  RESULTS  AND  DISCUSSION 

An  expansion  shockwave  tube  was  used  for  all  of  the  experiments 
cited  in  this  document.  The  booster  section  had  an  interior  dia¬ 
meter  of  80  mm,  it  was  2000  mm  long;  the  passive  booster  stage  had  an 
interior  diameter  of  185  mm,  and  it  was  5000  mm  long.  They  both 
used  the  same  type  of  shockwave  tube  diaphragm.  The  operational 
fuel  used  for  the  booster  section  and  the  passive  booster  was  ambient 
nitrogen.  Attached  to  the  nozzle  port  was  a  1.22  m  slotted  piezo¬ 
electric  sensing  device,  used  to  record  the  pressure  over  time  curve 
data,  as  outlined  in  Figure  4.  From  an  analysis  of  the  curve  data 
in  the  fiture,  we  can  make  the  following  assumptions: 


(1)  Figures  4a  and  4b  are  split  into  M  =  2.29  and  M  =2.22 

s  s 

for  the  pressure  curves.  Since  the  forward  platform  is  lower  than 
the  lower  platform  (p.  <  p,,  ) ,  we  can  diagnose  this  as  pre-expansion 

Z  JD 

flow. 

(2)  Figure  4c  is  for  the  pressure  cufves  at  Mg  =  2.11.  The 
forward  and  aft  platforms  are  nearly  equal  (  P,vp^  ) ,  so  we  can 
recognize  this  as  a  condition  of  near  "greatest  disturbance"  for 
the  compression.  The  following  data  is  revealed  from  the  convex 
curvature  that  manifests:  under  those  conditions  of  near  "greatest 
disturbance"  in  compression,  the  upper  reaches  of  the  boundary  layer 
did  not  smooth  out  as  theoretical  predictions  had  expected.  Further 
more,  there  was  a  turbulence  present  which  can  only  lead  us  to 
conclude  that  it  was  a  result  of  mutual  interference. 

(3)  Figures  4d,  e  and  f  are  split  into  curves  for  M  =  1.94, 

b 

M  =  1.84  and  M  =1.73.  The  forward  platform  was  higher  than  the 
s  s 

aft  platform  (p2  <  indicating  a  post  compression  flow.  The 

time  required  to  pass  through  both  platforms  was  about  3  ms.  The 

2 

velocity  of  the  primary  shockwave  was  approximately  6x10  m/s.  The 

primary  shockwave  required  approximately  2-3  ms  to  reach  its  peak, 

and  the  secondary  shockwave  required  about  2-4  ms.  The  velocity 
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of  the  secondary  shockwave  was  1°  t's.  I£  we  consider  the  seconder 

shockwave  as  a  primary  wave  of  15.5  r.tn,  the  time  required  for  it  to  pass 
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through  the  sensing  devices  would  be  aocuc  iO  ms;  the  primary 
shockwave  passed  through  the  boundary  conditions  and  a  division  was 
formed  that  is  demarcated  with  the  wide-line  recordings.  The  time 
required  for  the  primary  shockwave  to  pass  through  the  sensing 
devices  was  approximately  10  ^  ms.  For  all  three  of  the  above  con¬ 
ditions,  the  elapsed  times  were  about  3  ms.  From  this,  we  can 
assume  the  secondary  shockwave  to  be  a  compression  wave. 

(4)  Figure  4g  is  for  the  curve  of  Mg  =  1.63.  The  aft  platform 
for  this  condition  was  not  clearly  measured,  so  the  pressure 
booster  stage  reaction  was  inconclusive. 


The  results  of  the  above  observations  and  analyses  can  be 
summarized  as  follows:  >  2.11  is  a  condition  of  pre-expansion. 

The  boundary  layer  itwelf  is  actually  a  part  of  the  rarefaction 
wave.  Mg  <  2.11  is  a  condition  of  post-expansion.  The  boundary 
layer  here  exists  as  a  part  of  the  secondary  shockwave,  as  was  pre¬ 
dicted  in  bibliographical  reference  [3] .  Actual  observations 
support  this  assertion. 

The  expansion  shockwave  tube  flow  fields  exist  in  four  differ¬ 
ent  forms.  In  analysis  and  calculation,  we  must  recognize  which 
of  these  conditions  is  affecting  the  flow  under  study,  and  then 
make  the  proper  choice  as  to  which  calculation  methodology  is  best 
suited  for  use.  Figure  5  depicts  the  realms  of  each  of  the  four 
types  (limited  to  our  study,  aside  from  the  calculation  formulas) . 
From  the  figure,  we  can  see:  post-expansion  flow  occupies  a 
specified  width  within  its  realm.  Many  of  the  actual  parameters 
used  are  derived  from  this  region.  If  post-expansion  exists  in 
the  secondary  shockwave,  then  its  calculation  method  will  be  diff¬ 
erent  from  that  for  pre-expansion.  If  the  post-expansion  exists  as 
a  compression  wave,  then  the  rarefaction  wave  is  produced  in  the 
ore-expansion  area  in  abnormal  fashion.  Parameters  used  for  fore 
and  aft  wave  systems  are  similar.  Because  the  post-expansion  area 
can  be  substituted  in  the  formula  for  pre-expansion,  the  resultsd 
are  found  in  a  very  convenient  manner. 

I  would  like  to  thank  Comrades  Gu  Shengxue  and  Li  Zhenhua  for 
their  help  in  taking  measurements  during  this  experiment. 
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OBSERVATION  OF  WAVE  DIAGRAMS  FOR  SHOCK 
TUBE  WITH  THE  DIVERGENT  NOZZLE  AT 
DIAPHRAGM  SECTION 

* 

Yu  Hongru,  Lin  Jianmin,  Li  Zhongfa,  Gu  Jiahua 
(Institute  of  Mechanics,  Academia  tinicat 

Abstract 

Piezoelectric  pressure  data  were  obtained  in  i  185mm  shock  tube 
with  divergent  section  at  shock  Mach  number  range  1.  5  <  M,  <  2.  5  - 
It  is  shown  that  in  uaderexpansioa  flow  exist  rarefaction  wave  as  pre¬ 
vious  predict,  but  in  overexpansion  flow  exist  compressive  waves  in¬ 
stead  of  secondary  shock  wave  which  was  predicted  by  Roscizewcki 
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A  NUMERIC  CALCULATION  METHODOLOGY  FOR  THE  INFLUENCES 
OF  GROUND  TO  AIR  INTERFERENCE 

Chen  Zhengcai 

(Chinese  Aerodynamic  Research  and  Development  Center) 

1 .  FOREWORD 

Take-off  and  landings  are  important  stages  of  flight  for  an 
aircraft.  When  an  aircraft  is  flying  very  close  to  the  ground,  it 
can  be  influenced  by  the  effects  of  ground  to  air  turbulence,  and 
this  is  called  the  "ground  effect".  Generally  speaking,  it  causes 
the  aircraft  lift  to  increase.  Normally,  the  angle  of  attack  for 
take-off  and  landings  is  a  =  8-12°,  and  for  some  foreign  aircraft, 
it  may  be  as  high  as  approximately  20°.  With  wing  flap  deflection 
angle  5  =  25°-45°,  the  ground  effect  can  be  quite  serious,  even  to 
the  point  of  producing  obviously  non-linear  effects.  The  longitu¬ 
dinal  moment  of  force  and  the  trim  can  both  produce  very  obvious 
changes.  As  a  result,  when  designing  new  aircraft,  the  topic  of 
ground  effect  has  become  an  absolute  necessity  for  inclusion  in  the 
research. 

To  determine  the  unique  ground-atmospheric  influences  on  new 
aircraft,  what  is  most  often  used  is  a  fabricated  experimental 
floor  model  for  use  in  a  low  velocity  wind  tunnel.  A  level  flow 
stream  across  the  floor  is  used  to  simulate  the  ground.  The  floor¬ 
boards  are  usually  permanently  attached  in  the  wind  tunnel.  Clearly, 
this  cannot  actually  simulate  the  true  conditions  because  in  true 
flight  conditions,  the  aircraft  is  what  is  really  moving.  The  ground 
does  not  move  relative  to  the  atmosphere.  Nevertheless,  in  testing, 
the  floor  plates  and  the  atmosphere  are  what  is  moved.  This  reveals 
the  reasoning  for  an  attached  surface  layer  existing  on  the  floor. 
This  can  be  overcome  by  using  the  atmospheric  wind  tunnel  flow  velo¬ 
city  to  account  for  floor  motion.  However,  because  of  the  structure 
and  other  factors,  it  is  not  really  adequate. 
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After  using  lift  linearization  and  precision  lift  plane 
(sectional)  theories  to  calculate  the  wing's  ground  effect,  the 
results  proved  to  have  relatively  large  error  differences,  even  so 
much  so  as  to  possibly  directly  oppose  the  experimental  outcomes. 
Later,  after  pouring  over  the  calculation  methodology  for  the  ground 
effect  in  lift  plane  theory,  we  saw  that  there  was  no  calculation 
results  applicable  to  the  aircraft  as  a  whole.  Therefore,  in 
research  on  the  problems  of  ground  effect,  this  has  generated  a 
great  deal  of  concern. 

In  the  studies  of  aerodynamic  calculations,  bibliographical 
reference  [1]  provides  an  introduction  to  the  methodlogy,  and  from 
it  we  can  extend  these  learnings  to  a  vast  range  of  applications, 
such  as  to  the  validity  of  complicated  calculations  or  even  to 
simpler  methodologies.  This  document  uses  some  of  those  methods 
to  calculate  the  ground  effect  for  the  wing  and  the  entire  aircraft. 
In  calculation,  the  planar  surface  and  tail  profile  are  both  con¬ 
sidered  as  being  non-linear,  the  flaps  and  the  fuselage  are  consi¬ 
dered  deflectors.  This  allows  the  results  of  calculation  to  be 
sufficiently  accurate. 

2.  CALCULATION  METHODOLOGY 

For  the  problem  of  ground  effect,  a  solution  for  air  flow 
pressure  can  be  found  in  the  equation: 

dx •  dyl  dj*  (1) 

In  the  equation,  x,  y  and  z  are  the  three  coordinate  axes  within 
the  coordinate  system.  j>  is  the  interference  velocity  potential. 

To  satisfy  all  of  the  boundary  conditions,  we  include  the  following: 
influx  flow  non-interference  and  wing  surface's  rear  edge  Kutta  con¬ 
ditions  as  well  as  the  planar  surface  and  ground  interference  flow 
stream  conditions.  The  problems  associated  with  these  unique  con¬ 
ditions,  in  comparison  with  the  calculations  for  normal  flight 
aerodynamics,  are  that  we  cannot  bypass  the  ground  flow  conditions. 
Because  of  the  limitations  of  the  ground  surface,  we  handled  this 
as  an  infinitely  large  "mirrored  surface".  By  using  the  mirror 


surface  profile,  we  are  able  to  satisfy  the  special  conditions  for 
this  problem. 

This  document  separates  the  aircraft  calculations  into  segments 
for  the  wing  system  (wing  and  elevators)  and  the  rotation  around  the 
body's  planar  surface.  A  finite  element,  fundamental  solution  prob¬ 
lem  solving  method  is  used  for  the  equations.  At  the  same  time,  the 
fundamental  solution  figures  for  a  calculation  model  of  the  air¬ 
craft's  profile  are  given  below: 


Figure  1  Several  models  of  aircraft  wing  ground  effect 
calculations. 

1 — fundamental  solution  of  vortex?  2 — control  point;  3 — aircraft 
( illeg . ) ;  4 — ground;  5 — image  of  wing?  6 — (a)  linear  lift  model; 
7  —  fundamental  solution  of  vortex  image;  8— (b):  lift  model 
linear  boundary  conditions);  9 — (c) :  lift  model  (non-linear 
boundary  conditions) 


WING.  The  aircraft  wing  is,  along  its  level  surface  area, 
sectioned  off  into  trapezoidal  grid-.  ^i-tribu  ‘red.  a -e  =  .r  I .  mesh 
point  in  the  grid  is  a  fundamental  hurrc  solution  for  the  successive  1- 
shaped  vortex  intensities  and  the  consecutively  stronger  vortex 
intensities  which  are  similar.  The  image  of  the  fundamental  solution 
intensity  is  generally  the  same  as  the  wing  image's  fundamental 
solution.  Only  the  axis  of  the  vortex  will  differ  (see  Figure  1). 

The  boundary  condition  control  point  is  established  from  the  center 
of  the  grid.  Thus,  the  problem  can  be  solved  in  the  matrix  equation 
below: 

C A-  Al{r)  =  -iV.)  (2) 


In  the  equation,  A  (A)  is  the  wing's  fundamental  solution  (maged 
fundamental  solution)  coefficient  of  influence.  T  is  the  funda¬ 
mental  solution  for  intensity.  is  the  planar  orientation 

velocity. 
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Figure  2.  Tail  vortex 
r — -*i mate  system  '  .rug 

*  --  g ;  b-“pcinL  „ . .  _~  rur  ectcr^ 
of  wing;  c — tail  vortex;  d — ground 
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Figure  3.  Quantitative  lift 
for  two-dimensional  wing 
ground  effect, 
a — non-linear  conditions; 
b — linear  conditions 
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Figure  5.  Two-dimensional  wing 
ground  effect  lift  coefficents 
which  follow  changes  in  ground  to 
air  distances. 

a — linear  boundary  conditions; 
b — non-linear  boundary  conditions 
c — free  form  accuracy  figures 
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Figure  4.  Two-dimensional 
wing  with  and  without  We  are  obliged  to  point  out  the 

ground  interference  load  ,  , ,  . ,  . ..  ..  , 

distributions.  following;  the  problem  with  the  ground 

a — free  form;  b — linear  effect  calculations  using  linear  bound- 

c--non-l inear 'boundary  ary  conditions  is  that  the  calculation 

conditions  error  can  be  fairly  large.  This  docu¬ 

ment  uses  non-linear  boundary  conditions  (see  Figure  1).  At  the 
same  time,  the  U-shaped  vortex  of  the  tail  eddy  at  the  rear  edges 
of  the  wing  are  dragged  through  the  boundary.  Furthermore,  for 
the  unlimited  planar  surface  at  the  aft  section,  there  is  a  gradual 
linearization  to  ground.  An  equation  for  the  principles  of  the 
tail  vortex  (Figure  2)  is; 


BODY.  Suppose  within  a  body  of  revolution,  the  axis  line  is 
divided  into  a  dispersion  source  and  a  vortex  origin,  wherein  the 
fundamental  solution  can  substitute  for  the  interference  field. 

The  source  lies  within  the  grid's  center;  the  vortex  origin  lies 
at  the  head  of  the  grid.  Opposite  the  ground,  there  is  no  change 
in  the  source  image  and  the  intensity  of  the  vortex  origin  is 
equivalent,  but  it  lies  opposite  the  axis  line.  The  boundary  cond¬ 
ition  points  are  taken  in  a  source  location  that  is  opposite  the 
planar  surface,  either  upper  or  lower  profile.  The  matrix  solution 
can  be  found  in  the  equation: 

ZB-b-G- Gl‘0-R}=  (4) 

In  the  equation,  3  and  G  express  dispersion  and  vortex  source 
influence  coefficients.  Q  and  R  differentiate  between  their  res¬ 
pective  intensities.  "1"  expresses  the  image  quantity. 

3.  CALCULATION  RESULTS 

Calculation  results  for  the  aircraft  components  and  the  air¬ 
craft  as  a  whole,  as  mentioned  above,  are  given  below. 

TWO  DIMENSIONAL  WINGS 

This  document  uses  an  aspect  ratio  of  X  =  80  for  the  rectangular 
wing  and  the  chord  orientation,  which  is  separated  into  two  sections 
for  two  dimensional  results.  This  seems  to  be  the  most  rational 
methodology. 

Ground  effect  generally  causes  conventional  changes  that 
increase  aircraft  wing  lift.  The  reasons  for  this  are  shown  in 
Figure  1(a)  which  illustrates  the  linear  lift  theory.  However,  the 
fact  of  the  matter  is  that  there  are  times  when  lift  may  actually 
be  reduced  (see  Figure  3).  Figure  1(b)  illustrates  a  model  that 
considers  the  chord  angle  to  load  distribution.  The  rationale 
here  is  that  of  Figure  1(c)  which  illustrates  a  model  with  non-linear 
boundary  conditions.  This  is  the  calculation  methodology  that  is 
used  in  this  document.  c . 


Figure  4  shows  chord  orientation  to  load  distribution,  and 
Figure  5  illustrates  when  a  -  4°  lift  coefficient,  which  follows 
the  curvilinear  changes  in  distance  and  altitude  h.  As  one  can 
see,  the  linear  boundary  conditions  and  non-linear  boundary  condi¬ 
tions  have  large  differences  in  their  calculation  results.  When 
the  distance  between  ground  and  altitude  is  equal  to  one  times  the 
chord  length  value,  the  ground  effect  has  already  become  very 
feeble.  Lift  is  increased  approximately  5%  and  at  a  point  which  is 
two  times  the  chord  length  value,  it  is  minimal  and  may  even  be 
neglected  since  it  is  not  worth  the  trouble  of  calculating.  However, 
when  it  is  lowered  to  0.5  times  the  chord  length,  it  is  very  clearly 
seen  and  lift  shows  a  quantitative  increase  of  as  much  as  15%. 

THREE  DIMENSIONAL  WINGS 

Figure  6  and  Figure  7  are  for  X  =2.8  for  a  rectangular  wing's 
lift  and  moment  of  force  calculation  results  [2],  At  the  same  time, 
they  provide  test  values  and  results  from  international  experiments 
which  can  be  used  to  make  quantitative  comparisons  with  our  results 
to  adequately  satisfy  our  needs.  Moreover,  the  ground  effect 
causes  the  linear  lift  slope  to  increase  =  0.2c  (c  is  the  average  * 
aerodynamic  chord  length)  when  it  is  greater  than  90%.  When  h  = 

1.0c,  within  10%  approximately,  we  can  see  the  three  dimensional 
wing  ground  effect  is  a  much  more  severe  problem  than  when  compared  - 
to  ground  effect  for  a  two  dimensional  wing.  This  is  because  of 
the  greater  image  interference  of  the  tail  vortices  system. 

ENTIRE  AIRCRAFT 

The  calculations  for  the  entire  aircraft  were  done  as  outlined 
in  bibliographical  reference  131  of  this  document.  The  model  for 
the  body  was  a  12  scale  replica  tapered  body.  The  wing  was  tri¬ 
angular,  X  =  3.0.  The  tail  was  waist  tapered,  X  =  4.0  (Figure  8). 

Figure  9  depicts  a  curve  of  the  entire  aircraft's  moment  cf 
force  characteristics.  The  test  values  of  the  figure  were  pro¬ 
vided  from  bibliographical  reference  [3]  which  were  taken  from  the 
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ground  interference  lift 
characteristics. 


a — point  used  for  moment  of 
force;  b — test  figures;  c— 
calculation  figures  of  biblio¬ 
graphical  reference  [2] ;  d — 
calculation  figures  of  this 
document 
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Figure  7.  Rectangular  wing 
(\=2.8)  with  and  without 
ground  interference  moment 
of  force  characteristics, 
a— test  figures;  b — calculation 
figures  of  bibliographical 
reference  [2] ;  c — calculation 
figures  from  this  document 
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Figure  8.  Calculation  mcual 
for  the  entire  aircraft, 
a — consideration  point  for 
moment  of  force;  b — center 
line;  c — ground 


c., 


Figure  9.  Entire  aircrai  :  vitr. 
and  without  ground  interference 
moment  of  force  characteristics, 
a — test  figures  (3);  b--calculation 
figures  of  this  document 
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wind  tunnel  floor  test  results.  The 
-ail  was  stabilized  with  a  deflection 
of  0.2°,  which  was  the  same  as  for  the 
other  tests.  Regardless  of  the  pre¬ 
sence  of  ground  interference,  the  cal¬ 
culation  values  and  the  experimental 
values  both  coincided  quite  nicely. 

Figure  10,  when  the  wing  flap 
deflection  is  6  *  40°,  shows  moment 
of  force  curves  for  the  whole  aircraft. 
Although  there  is  a  large  flap  deflec¬ 
tion,  it  is  due  to  the  non-linear 
boundary  conditions  which  are  used  in 
this  document,  and  the  calculated 
figures  coincide  closely  to  other  [3] 
test  results. 

The  interference  between  the 
ground  surface  and  the  aircraft  obvious¬ 
ly  changes  with  the  moment  of  force 
curvilinear  slope,  and  it  influences  the  trim  characteristics  as 
well . 

4.  MAJOR  CONCLUSIONS 

1.  For  aircraft  very  close  to  the  ground,  when  using  a  linear 
boundary  condition  to  calculate  all  of  the  results  for  ground 
effect,  the  error  difference  can  be  very  large.  We  must  utilize 
non-linear  boundary  conditions  that  will  allow  us  to  account  for 
conditions  of  deflection. 

2.  Analysis  in  basic  linear  lift  theory  allows  us  to  obtain 
ground  effect  which,  in  general,  causes  aircraft  wing  lift  to 
increase.  Because  the  flaps  have  specific  deflections  at  different 
altitudes,  ground  interference  can  actually  cause  a  decrease  in 


Figure  10.  Entire  aircraft 
(5=40°)  with  and  without 
ground  interference  moment 
of  force  characteristics, 
a — test  figures;  b — calcul¬ 
ation  results  of  this 
document 


lift.  Linear  lift  theory  does  not  account  for  these  conditions. 


3.  The  ground  to  aircraft  flight  aer-odynamics  are  severely 
influenced  by  ground  turbulence,  and  ground  turbulence  clearly 
alters  the  lift-line  slope  and  in-flight  aerodynamics. 

4.  The  finite  element  fundamental  solution  calculation  method¬ 
ology  is  used  to  calculate  flight  ground  effect  and  is  a  very 
effective  method. 
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A  CALCULATION  METHOD  OF  GROUND  EFFECTS 
FOR  THE  AIRCRAFT  ' 

Chen  2hengcai 

(China  Aerodynamic  Research  and  Development  Centre) 

Abstract 

When  an  aircraft  js  flying  near  the  ground,  it’s  aerodynamic  charac¬ 
teristics  may  be  influenced  by  the  ground.  Using  the  finite  element 
solution  method  and  the  linear  and  nolinear  boundary  conditions,  this 
paper  gives  a  calculation  method  for  the  ground  effects.*  Here  the 
ground  is  considered  as  an  infinite  image  plane  and  both  the  flap  and 
the  elevator  can  deflect  too. 

The  calculation  results  for  two-or  three-dimensional  wing  and  the 
wing-body-tail  combinations  are  presented  and  these  results  are  in  a 
good  agreement  with  the  experiment  data  or  other  theoretical  results. 

The  calculation  results  for  deflecting  control  surfaces  are  only  given 
by  our  method. 

The  calculation  results  proved  that  using  the  linear  boundary  condition 
and  the  lifting  line  theory  may  sometimes  produce  the  wrong  results, 
especially  as  the  wing  is  guite  near  the  ground.  By  means  of  the 
aolia  ear  boundary  condition  and  the  lifting  surface  theory  this  method 
can  get  good  results. 
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ON  DETACHED  SHOCK  NAVES  OF  SPHERICAL  BODIES  AT 
TRANSONIC  VELOCITIES 

Zheng  Zhichu 

(Institute  of  Mechanics,  Academia  Sinica) 

1 .  FOREWORD 

In  recent  times,  the  primary  orientation  of  research  for  thick 
bodied  transonic  linear  flow  has  been  within  the  realm  of  :  . 

Within  the  range  of  0.so<M,<i.05  ,  especially  for  the  area  closest 

to  i  ,  results  of  theoretical  calculation,  experimentation  and 

wind  tunnel  tests  for  spherical  bodies  have  provided  evidence  that 
there  are  a  great  many  differences  which  exist  across  a  wide  range 
of  detached  spherical  bodies.  For  an  accurate  measurement  of  the 
shock  waves  and  certain  changes  they  undergo  as  the  Mach  number 
changes,  there  is  still  no  prescribed  best  method. 

The  calculations  referred  to  in  bibliographical  reference  [1] 
point  out  the  following:  within  the  vicinity  of  M.s:,  the  changes 
of  the  detached  shock  wave  as  it  follows  the  changes  in  the  Mm 
values  are  even  more  sensitive  when  compared  to  the  physical  plane 
surface  pressure  distributions.  As  a  result  of  this,  at  transonic 
velocities,  it  is  more  accurate  to  use  the  spherical  body  detached 
shock  wave  measurements  as  the  Mach  number  changes  instead  of 
using  the  pressure  coordinates. 

Based  on  the  methods  outlined  in  bibliographical  reference  [2] , 
we  need  only  follow  the  rules  outlined  for  detached  shock  waves 
which  follow  value  changes  within  the  vicinity  of  .  We 

can  use  a  small  arms  ballistic  range  and  common  precision  instru¬ 
ments  for  measurement  to  obtain  a  resistance  coefficient  for  the 

This  document  was  received  in  1983,  revisions  were  received  on 
February  29,  1984. 


spherical  body  at  the  vicinity  of  :  .  Afterwards,  following  a 

comparison  with  generalized  wind  tunnel  results,  free-f light  data 
or  magnetic  suspension  modeled  wind  tunnel  results,  we  can  research 
the  influence  of  the  tunnel  wall  and  center  of  resistance  in  regard 
to  spherical  body  aerodynamic  coefficients  [3] . 

In  consideration  of  the  aforementioned,  it  is  of  great  signi¬ 
ficance  to  note  that  all  of  the  data  provided  within  the  realm  of 
transonic  velocities  when  compared  to  the  data  for  the  region  of 
'  '  :  ,  runs  counter  to  the  simple  equation  for  detached  shock  waves 
of  spherical  bodies. 


Since  there  is  no  center  of  resistance  interference  present  in 
a  ballistics  test  range,  it  is  comparable  to  a  modern  transonic 
velocity  wind  tunnel.  The  ballistics  test  results  are  quite  similar 
to  flight  results  in  comparison  14] .  This  document  uses  the  results 
of  ballistics  range  experiments  as  a  foundation,  then  a  comprehen¬ 
sive  analysis  is  done  to  determine  the  influence  of  each  variable 
on  shock  wave  form  with,  for  instance,  the  rules  for  the  detached 
distance  6,  the  shock  wave  slope  tgf?,  and  the  shock  wave  post-sonic 
velocity  position  9*,  etc.,  which  follow  the  changes  in  the  Mach 
number  M  .  Furthermore,  with  all  of  the  assimilated  data  results  for 
experiments  conducted  in  recent  history,  using  those  values  in  the 
range  of  ,  a  comparison  of  actual  flow  for  spherical 

detached  shock  waves  of  various  forms  can  be  used  to  create  a  simple 
equation. 

2.  DETACHED  SHOCK  WAVE  FORMULA 

Spheres  take  the  form  of  symmetrical  axis  bodies.  Research 
of  spherical  detached  shock  wave  forms  only  requires  consideration 
of  the  meridian  plane  (see  Figure  1).  We  all  know  that:  when 
v, ,  the  physical  body  of  the  circular  arc  of  the  detached  shock 
wave  is  stationed  directly  in  front  of  the  sphere?  moreover,  in 
the  range  i-M.  when  it  is  within  the  region  of  -*■  1,  the 

spherical  detached  shock  wave  should  follow  the  values  according 

60 


Figure  1.  Polar  coor¬ 
dinate  system  for 
spherical  detached 
shock  waves. 


to  specifically  defined  laws.  A  single 
value  can  be  obtained  through  a  linear 
intersection  of  the  arc.  We  hypothesize 
the  use  of  a  curve  ratio  from  the  center 
point  of  the  sphere,  sectioned  in  conical 
fashion  to  express  shock  wave  form.  Its 
polar  coordinates  are  expressed  in  the 
following  equation  as: 

r  _  P _ 

~R  i-ecoss 


wnen 


sn  8  «  0°, 
T  =  _P 
R  I  -  £ 


and  by  using  substitution  in  the  above  equation,  we  get 
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Differential  equation  (1)  can  be  used  to  get  shock  wave  inclin- 
ination  B  or  the  corresponding  relationship  of  polar  coordinate 
angle  8 

••!>/-'  .  ,  „  J  (  i  ) 


Using  the  slope  of  the  shock  wave  to  create  a  relationship  of  post- 
M  shock  wave  values  where  M-post  =  1,  we  can  obtain  post  sonic 
shock  wave  polar  coordinate  angle  9*  and  e,  and  their  relationship 
is  as  expressed  below: 


so  we  have 


s ;  r. 7 *  =  s  cosd* '  s:r.^*  -  '  si n'3*  -  ( 1  - 

i  \  £■ 


iff-  .1-v—  ]} 


From  equations  (1;,  (3)  and  (4),  we  can  see:  rQ/R  £  are 

both  following  the  quantitative  changes  of  the  influx  Mach  number 

values  M  . 
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3.  AN. ANALYSIS  OF  THEORY  AND  EXPERIMENTAL  RESULTS 

For  an  analysis,  we  are  first  concerned  with  the  calculated 
values  for  5/R  and  9*.  We  start  out  with  the  experimental  results 
and  several  test  formulas,  and  we  selected  5  *  i/v  M1_- !  to  make  a 
self  changing  variable.  By  tracing  the  variables  5/R  and  9*,  we 
see  that  they  follow  the  changes  along  the  curves  ? ,  and  they  can 
be  differentiated  in  Figures  2  and  3.  From  the  figures,  we  can  see 
after  using  £  to  create  an  abscissa,  in  the  time  of  Mw  -  1,  we  can 
clearly  observe  the  current  variational  results  for  all  trends  of 
5/R  and  9*  as  they  follow  the  numeric  changes  of  the  M^.  In 
Figure  2,  as  an  example,  the  results  for  5/R  seem  to  appear  as 
"one  objective  to  find,  three  paths  to  take".  The  "one  objective" 
is  for  the  various  results  after  M  >  1.5  to  fall  in  place.  It 
seems  clear  that  in  hypersonic  velocity  ranges,  theory  and  test 
results  coincide  in  excellent  fashion.  The  "three  paths"  are  for 
the  numeric  values  calculated  after  M  <1.3,  the  test  results,  and 
the  test  formulas,  in  all  three  groupings  which  have  relatively 
great  differences  across  their  range  of  results. 

The  first  group  of  results  (5)  are  taken  as  wind  tunnel 
experimental  value  substitutions,  and  the  test  equation  can  be 
expressed  as: 

=  (5) 

A  O 

These  results,  at  M  >1.2  [8]  and  a  few  numeric  calculation 
results  seem  to  coincide  favorably  well.  At  M  =  1.05,  the  results 
of  other  calculations  are  similar. 

The  second  group  of  results  substituted  for  the  wind  tunnel 
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Figure  2.  Spherical  detached  shock  wave  distance 
following  numeric  changes  in  values  for  Mm . 
a — experimental;  b — numeric  solution;  c,d — sphere 
^.--formula  of  this  document 
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Figure  3.  Post-sonic  shock  wave  positions  following 
changes  in  the  values  for 

a,b,d,f,g,h,i,j, k — sphere;  c , e — spherocylindrical; 

1 — results  of  this  document;  m — numeric  solution 


(1),  (7),  also  have  a  simple  test  equation  that  can  be  used: 

-Jr  -0.156;  U  -1  (6) 

A 

In  the  range  M»  =  1  .-3—1  •  ,  the  aforementioned  results  from  (1), 

(7),  (8)  and  (9)  are  put  in  the  numeric  calculation  results  and 

they  compare  very  closely. 

The  third  group,  for  (10)  and  (4) ,  are  ballistics  range  results 
and  substituted  wind  tunnel  results.  When  Mm  >  1.05,  they  compare 
very  closely  to  the  results  obtained  from  the  firsg  group.  When 
<  1.05,  the  outcome  falls  somewhere  between  the  results  from 
the  first  and  second  groups.  Based  on  the  fact  that  the  ballistics 
range  equipment  is  clearly  used  only  as  a  comparison  flow  field  in 
lieu  of  actual  flow  and  its  unique  characteristics,  we  recognize 
that  the  third  group  of  results  run  counter  to  the  true  trends  of 
change  for  6/R  which  follow  the  variations  in  M  .  Furthermore, 
the  reason  for  creation  of  the  third  group  of  results  is  as  follows: 
within  the  vicinity  of  M_si  ,  calculation  of  numeric  values  and 
experimental  results  of  wind  tunnel  measurements  are  both  very 
difficult  to  accurately  ascertain.  For  the  many  various  types  of 
mistakes  that  could  be  made  in  obtaining  the  results,  a  self 
correcting  evaluation  result  was  needed  to  prove  accuracy.  The 
utilization  of  the  Mach  number  counter  flow-field  results  were 
assimilated  for  a  composite  comparison. 


As  a  foundation  for  the  third  group  of  results,  we  chose  the 
following:  ,=  . 

i '  f  ■ 

"  -  (7) 

=  0 . 4d£  '  (  8 ) 


By  using  formulas  (1),  (2),  (7)  and  (8),  we  can  calculate  <5/R,  3, 

r/r  ,  and  r/R  for  a  comprehensive  examination  of  the  results  com¬ 
pared  with  the  variations  in  the  Mach  number  values.  The  comparison 
conditions  for  the  experimental  results  and  experimental  formulas 
may  not  actually  reflect  the  conditions  for  true  flow. 


4.  RESULTS  AND  DISCUSSION 


Figures  (2),  (3),  (4),  (5)  and  (6)  differentiate  between  the 


various  value  changes  of  6/R,  9*,  r/rQ  and  r/R  as  they  follow 


relational  changes  in  the  values.  From  Figure  (2),  we  can  see: 
the  5/R  value  calculated  using  (7),  (3)  in  the  range  of  1  <  <4  for  the  thi 

group  of  results  compares  very  well.  In  comparison  with  the 
various  methods  currently  used,  in  the  vicinity  of  V,r  ; ,  the 
results  here  in  this  document  actually  are  at  odds  with  the  rules 
set  for  the  detached  shock  wave  distance  which  follows  changes  in 


the  values. 
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From  Figure  (3),  we  can  see:  the  results  of  9*  within  the 
ranges  of  1  <  <  4,  obgained  in  formulas  (7)  and  (4)  in  the 

third  group  of  equations  compare  very  well.  They  also  point  out  a 
few  important  manifestations.  Foremost,  at  Nt,  i ;  vicinity,  the 
value  of  9*  is  very  sensitive.  The  differences  in  variational  e 
can  greatly  influence  the  results  of  calculation.  For  example,  if 
we  first  use  formula  (7)  and  the  variational  z  as  expressed  in  the 
example  formula  below:  } 

wU-.-'iY-  •’  (9) 

we  can  start  our  comparison.  From  Figure  (3),  one  can  see  that 
the  obtained  value  for  6  *  in  formulas  (7)  and  (8),  within  the 
range  of  Mm  >  1.2,  is  pretty  much  the  same.  However,  at  <  1.05, 
there  is  a  great  discrepancy.  The  variational  e  obtained  in 
formula  (7)  showed  9*  results  from  the  third  group.  Furthermore, 
with  the  variation  from  formula  (9),  at  M  =  1.01,  9*  cannot  be 
solved  at  all.  Tbis  occurrence  could  explain  why,  at  conditions 
of  hypersonic  velocity,  the  calculation  methods  show  a  wide  range 
of  difference  in  their  results,  especially  at  <  1.05,  where  it 
seems  to  be  most  difficult  to  figure.  This  is  the  reason  the  shock 
wave  form  does  not  accurately  follow  the  rules  for  changes  in  the 


M  values. 
00 


Additionally,  the  formulas  in  this  document  obtained  values 
for  9*.  At  M  -  1,  there  is  a  gradual  decrease,  and  at  M  = 


V,*? 


1.000001,  9*  54°.  It  along  with  formula  (11)  linear  flow  stream 
at  =  1  range,  accounts  for  ever,  linear  flow  of  the  sphere  at 
sonic  speeds.  From  the  values  for  M_— ;  ~  and  ,  we  can 

solve  for  the  linear  course  at  sonic  velocity.  When  ■  :  ,  the 

detached  body  of  the  shock  wave  stabilizes  to  that  of  the  shock 
wave,  and  the  flow  following  the  wave  is  subsonic.  When  the  shock 
wave  is  in  extension,  it  is  limited  by  the  value  of  the  Mach  number, 
and  the  flow  around  it  is  hypersonic.  As  a  result,  ohe  shock  wave 
must  have  a  sonic  point  along  its  curve,  as  well  as  a  value  for  9* 

<  90°  when  values  M  -*■  1  are  present.  Also,  when  M  -  1,  there  is 
no  limit  to  the  distance  of  the  physical  plane  of  the  shock  wave. 
This  is  because  as  M._  -►  1,  the  sonic  points  along  ^hock  wave 

are  in  a  forward  position. 


*  □  M  ■  1  .01 

—  *  !•??- . 

,  ■  M  - 

A  M  ■  1  .*5 


o  :a  w  9« 

0  i0  «0  3? 

Figure  4.  Comparison  of  results  for  r/r  fell" wing 
the  changes  for  9 . 
a — formula  of  this  document. 

From  the  calculations  of  r/rQ  obtained  in  formulas  (2)  and 
(7),  Figure  4  shows  experimental  results  and  numeric  figures  that 
coincide  in  excellent  fashion.  We  should  point  out  that  the  rQ/R 
values  are  counter  to  the  results  for  variation  e.  Figures  (5) 
and  (6)  show  that  within  the  range  of  1  <  Mw  <  4,  results  of  r/R 
follow  the  changing  values  of  Mm.  From  the  figures  above,  we  can 
see  that  the  numeric  solution  for  it  and  the  experimental  results 
coincide  quite  nicely. 
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Figure  5.  Detached  shock 
wave  form  of  spherical 
body  at  transonic  velocity, 
a — formula  of  this  document 


Figure  6.  Form  of  detc-o  'ar  s..ook 
wave  of  sphere  as  it  follows  che 
changes  of  values  for  Ma. 
a — formula  of  this  document 


In  summation  of  the  aforementioned,  this  document  provided  a 
formula  for  calculations  within  the  range  of  1  <  <  4.  Within 

the  vicinity  of  =  1,  actual  comparisons  are  at  odds  with  the 
changing  values  which  the  spherical  detached  shock  wave  should, 
by  rule,  follow.  The  results  can  be  applied  to  slotted  wall  inter¬ 
ference  in  transonic  wind  tunnels  and  in  observations  of  experimental 
booster  flow  fields.  It  is  also  applicable  to  numeric  calculations 
and  problem  solving  analyses  methodology  which  provide  results 
used  in  making  comparisons.  In  regard  to  spherocylindrical  and 
semi-conical  spheres  of  a  <  8°,  calculations  of  di-sonic  velocity 
detached  surface  shock  wave  bodies  and  experimental  results  coincide 
favorably  well  in  comparisons. 


The  author  would  like  to  express  his  sincere  appreciation  to 
Professor  Lin  Tongjia  for  his  guidance. 
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ON  DETACHED  SHOCK  WAVE  OF  SPHERE  MOVING 
WITH  TRANSONIC  VELOCITIES 


Zheng  Zh  ichu 
of  Mechanics,  Academia  Scienca- 


Abstract 

la  recent  years,  there  has  been  an  increased  interest  in  the  shock 
shapes  about  blunt  bodies  at  M.asl.  The  experimental  data  which  were 
obtained  by  testing  the  sphere  in  oallistic  range  facility  have  been  an¬ 
alyzed  and  a  signle  expression  for  the  shock  wave  profiles  of  sphere 
is  presented  in  this  paper.  The  polar  form  of  this  equation  is 

TT"i1  +  xXi*«5s-  wbere  ^‘“ensionless  value  of  shock  detach- 


®ent  distance  ^  and  variation  £  are  function  of  uniform  stream  Mach 


number  M„.  In  the  1<  M.  <4  range,  especially  at  M.asl,  dimension¬ 
less  quantities  — ^  ,  9*  and  — —  which  were  calculated  by  equation 


given  in  this  paper  show  excellent  agreement  with  free  flying  data. 
Also,  this  equation  can  be  extended  to  calculate  the  shock  wave  profi¬ 
les  of  hemisphere-cylinder  and  spherically  blunted  cones  (half  angle  a<8') 
Furthermore,  results  in  this  paper  can  be  extended  to  study  the  wall 
disturbance  and  test  section  flow  field  in  transonic  wind  tunnel. 


SEPT.  1984 


JOURNAL  OF  AERODYNAMIC  SCIENCE 


I 

i 

FINITE  DIFFERENCE  COMPUTATIONS  FOR  RADIAL  WING  FLOW  IN  A 
|  TRANSONIC  TWO  DIMENSIONAL  SLOTTED  WALL  WIND  TUNNEL 

Zhang  Naiping 

(Northwestern  Polytechnical  University) 

1 .  FOREWORD 


Experimental  data  in  a  wind  tunnel  is  influenced  by  a  great 
many  factors,  chief  among  which  is  the  influence  of  slotted  wall 
interference.  This  is  a  very  serious  factor  in  a  transonic  velo¬ 
city  field.  By  using  a  di-slotted  or  quadra-slotted  transonic 
velocity  wind  tunnel  duct,  we  can  cause  a  large  reduction  of  this 
particular  type  of  interference.  However,  when  the  Mach  number  is 
relatively  high,  it  makes  little  difference  whether  in  theoretical 
research  or  experimentation  because  it  is  just  as  difficult  to  cope 
with  the  interference  in  either  case. 

The  classical  tunnel  wall  interference  theory  considers  the 
influence  of  compression  characteristics  in  linear  subsonic  velocity 
as  its  foundations.  However,  when  the  airfoil  surface  is  primarily 
in  hypersonic  velocity  fields,  a  linearized  theory  cannot  accurately  - 
provide  pressure  values.  Moreover,  the  shock  wave  displacement 
created  by  influences  of  the  tunnel  wall  cause  changes  to  occur  in 
the  pressure  coefficients.  Also,  a  linearized  theory  cannot 
account  for  these  changes  in  the  form  of  an  equation. 

Here  and  abroad,  we  have  already  applied  the  finite  difference 
method  and  relaxation  techniques  in  order  to  find  solutions  for 
radial  wing  flow.  It  has  also  been  used  within  radial  wind 
tunnels  for  solutions  to  airfoil  flow  problems  11,2,3)  mainly  in 
cavity  wall  applications.  This  document  makes  use  of  these  methods 
to  calculate  the  upper  and  lower  slotting  bypass  flow  at  transonic 
velocity  for  radial  wings  in  a  bypass  wind  tunnel.  The  transonic 


wind  tunnels  in  the  Research  Institute  of  Aerodynamics  and  Gas- 

2 

dynamics  at  Stuttgart  University  (IAG)  260x700  mm  and  in  the 

Institute  of  Aerodynamics  at  Northwestern  Polytechnical  University, 
2 

100x300  mm  ,  were  chosen  as  the  two  computational  examples.  Only 
blockage  interference  conditions  at  zero  angle  of  attack  are  cal¬ 
culated  herein. 

2.  EQUATIONS  OF  VELOCITY  POTENTIAL  AND  BOUNDARY  CONDITIONS 

The  calculations  of  this  document  use  the  applied  two-dimen¬ 
sional  transonic  velocity  small  perturbation  velocity  potential 
equation,  which  is 

<1  - -?>„  =  o  (1) 

This  equation  is  a  mixed  non-linear  equation.  After  application 
of  the  small  perturbation  conditions,  we  have 

and  herein,  M  is  the  primary  Mach  number. 

After  solving  equation  (1) ,  all  of  the  necessary  boundary 
conditions  for  the  wing  are  expressed  as: 

(3) 

In  the  equation,  F(x)  is  a  function  of  the  wing  exterior  surface; 
lower  coordinate,  u,  0  differentiates  among  expressions  of  lower 
and  upper  wing  surface.  Furthermore,  the  Kutta  boundary  conditions 
are 

9,{x,  +o>“f.(x, -0)  (4) 

9i(*t  +0)  —  0)  (5) 

In  the  equations,  +0  and  -0  are  differentiated  to  express  upper 
and  lower  placement  of  wing  trailing  edge  flow  expansion  velocity 
potential.  For  example,  if  H  replaced  the  trailing  edge,  then  we 
would  have 

(6) 

9<x,  +0)  —  p(.x,  -  0)  +0)  -?(**,  -0)  —  r 


This  document  was  received  on  Nov,  17,  1983,  revisions  were 
received  on  March  9,  1984. 


70 


Upper  and  lower  slot  conditions  in  the  wind  tunnel  can  be 
written  as 


c,cp 


c. 


d'<p 
dx  dy 


- 0 


(7) 


Equation  (7)  includes  different  types  of  boundary  conditions  in 
the  wind  tunnel.  In  regards  to  the  slotted  wall  adhesion  flow, 
we  havec.  =  o,c,*i,e,-i//?,  an^  the  above  equation  would  be  written 

as 


di T  1  os  i  o'x 

—  hi  ■  -  T  i 

dx  R  dy  dx  dy 


(8) 


;.-i|ncsc|'  -1  ± )  (9) 

n  V  2  d  > 

In  this  equation,  !  is  a  slot-width  t  and  interval  d  includes  its 
relevant  parameters.  In  regards  to  the  principles  of  a  slotted 
wall  wind  tunnel,  we  can  approximate  the  boundary  conditions  as: 


dip  ,  d'<p 

— -  T  ‘-Z - r—  =  J 

ox  dx  dy 


(10) 


The  above  equation  may  also  be  written  as 


(11) 


The  integral  equation  above  applies  upper  field  vernier  boundary 
conditions,  and  we  have 


(12) 

wherein  we  may  get  an  approximate  integral  form  for  the  boundary 
conditions. 


The  upper  and  lower  field  vernier  boundary  conditions  can  be 
reviewed  in  bibliographical  reference  II] .  If  the  wing  in  exper¬ 
imental  stages  is  placed  in  a  non-symmetrical  form,  slot  gap 
interval  can  be  differentiated  using  H^,  Hj  as  expressions.  When 
symmetrically  placed,  we  would  use  the  expression  =  Hj  =  H/2. 

We  also  have (p, Upper  and  lower  vernier  field  conditions 
are 

(13) 

lor  *(»+.) 
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Issue  3  Zhang  Naiping:  Finite  difference  computations  for 
airfoil  flow  in  two  dimensional  transonic  slotted 
_  wall  wind  tunnel 
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-e«(//,-y)-r  A*. 


upper  verruca! 


C* 

.  {*,  (H ,  t  ZH  i  —  y)  —  A?>, 

(' V .  .  (  3  „  .  lower  vortical 

r  “Vj  +  A#. 


(14) 


JA*.-*/.,-*,.,  (15) 

(Afl?;  =93,. 

2(//,  +  //,>' A*,  +  A*.)  (16) 

or  ?4  *  ini**'  ’*  **  ~  Vh- " -•  “  ) 

In  the  above  equations,  lower  coordinates  jt+0  and  Jt-0  replace 
upper  and  lower  vortical  surfaces. 


! 

3.  CALCULATION  PROCESS 

Equation  (1)  is  a  mixed  difference,  non-linear  form  partial 
differential  equation.  Within  the  field  of  flow,  we  must  first 
carry  out  a  velocity  discrimination  at  each  individual  mesh  point. 

In  regards  to  subsonic  and  sonic  velocity  points  <J>  and  <f>  use  a  - 
central  differences  equation.  For  those  points  at  hypersonic 
velocities,  we  use  rear-oriented  single  precision  differences. 

Aside  from  the  field  of  flow  boundary  conditions,  4>xx  is  also  taken  - 
from  among  central  differences.  With  the  line  relaxation  parallax 
y,  the  difference  equation  matrix  coefficients  take  on  a  tri¬ 
dimensional  linear  equation  form  which  we  can  use  as  a  tracking 
method  to  find  our  solution. 

The  slotted  wall  boundary  conditions  *  in  the  central 
difference  equation  embody  the  exterior  of  the  slotted  wall, 
namely  the  unspecified  points  (i,j)  and  (i,j+l).  From  the 
method  cited  below,  we  can  take  these  unspecified  points  and  obtain 
values  for  velocity  potential.  In  regard  to  the  upper  surface  of 
the  slotted  wall,  we  have 

■-  f.  (17) 

Ay* 
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1 


I  -  \  _  >  ' -i ' 

- r 


From  the  r.  ■-  .  value  in  these  :*c  equations,  we  get 

•-  -v'  •>-  /  (19) 

Herein,  for  the  lower  surface  of  the  slotted  wall,  we  have 

<r.  •  (20) 

with  the  slotted  wall  boundary  conditions  (12)  substituted  in  (19) 
or  (20) ,  we  have  „  , 

'> -377»"  (2D 


With  identical  wall  boundary  conditions,  for  the  wing  upper  and 
lower  surfaces,  we  can  obtain 


ty  2 

(<p„.  .  -J7T  (f-  1 


(<p„; Vi.  i»-«)  t  »  (?,)■.>*•« 


- - -  Ay,  - -  — Ay 

In  the  two  equations  above,  and  can  be  substituted  into 

equation  (3). 


Due  to  the  line  relaxation  axis  for  parallax  y,  each  mesh 
point's  difference  equation  can  only  emcompass  three  unknown  velocity 
potential  points,  so 

-o,  (25) 

In  the  equation,  coeff icients a,,b.,c,,d.  are  already  known  values. 


From  the  small  perturbation  hypothesis,  we  can  use  the 
following  for  the  calculation  of  the  pressure  coefficients: 

cy — fy-  (26) 

We  can  use  a  precision  formula  to  calculate  pressure  distributions 
along  the  wing  surface: 


C,m  kliL  ' 

-J 


cosa  y,  p,  $ma 


■*r~ 


1  >  (27) 


The  calculations  were  performed  using  the  Qitu  Jiate  University 
CDC  6600  computer. 
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4.  CALCULATION  RESULTS 


2 

The  calculations  for  the  IAG  260x700  mm  wind  tunnel  were  for 

a  non-symmetrical  model.  The  intervals  for  the  upper  and  lower 

2 

slots  m  the  model  were  separated  at  400  and  300  mm  .  Lift  was 
for  a  zero  angle  of  attack.  The  primary  purpose  of  this  document 
was  to  calculate  blockage  interference  at  zero  angle  of  attack. 


a  —  U'tkrmx’i'CXKtxM 


Figure  la.  Pressure  distribution 
:..;ves  far  a  NACA  0012  airfoil 
i-  differing  Mach  nos.  a--Langley 
"ssearch  Center  tesc  results;  In¬ 
finite  difference  calculations  of 
am s  paper 


Figure  lb.  Pressure  distri¬ 
bution  curves  for  a  RAE  104 
airfoil  at  differing  Mach 
numbers . 

a — finite  difference  calcula 
tions  of  this  paper? 
b — NPL  lab  results 


2 

A  RAE  1C4  airfoil  of  100  x  300  am  was  used  in  the  Northwestern  Polytechnical 
University  wind  tunnel.  A  NACA  0012  airfoil  was  used  in  the  IAG 
wind  tunnel.  The  calculation  results  of  this  document  were  com¬ 
pared  to  experimental  results  from  the  NPL  (National  Physical 
Laboratory)  in  England  and  from  Langley  Research  Center  in  the 
United  States. 

Figure  la  points  out  the  pressure  distribution  calculation 
results  of  Mach  0.73,  0.77  and  0.83  for  the  NACA-0012  airfoil. 

From  the  figure,  we  can  see  that,  aside  from  when  M  =  0.83,  other 
Mach  value  (0.58,  0.64,  0.68  results  were  not  included  in  the 
research  article)  results  and  those  from  the  Langley  Research 


Center  compare  favorably  well.  Clearly,  the  I AG  wind  tunnel  in 
the  range  of  these  Mach  values  has  the  least  blockage  interference. 
The  calculations  in  this  document  are  also  equally  reliable. 


Figure  lb  points  out  the  pressure  distribution  calculation 
results  at  Mach  0.6,  0.7,  0.85  and  0.87  for  a  RAE  104  airfoil. 
Moreover,  the  NPL  linear  flow  results  are  used  for  making  comparisons 
Aside  from  the  post  shock  wave  pressure  distributions  when  M  =  0.87, 
our  calculation  results  compared  very  closely  to  their  experimental 


results . 


Figure  2.  Coefficients 
of  biased  zero  lift  as 
produced  in  non- 
symmetrical  slotted  wall, 
a — slotting; 
b — experimental 


Figure  3.  Pressure  distribution 
curve  for  upper  and  lower  wind 
tunnel  slotting. 

a,c,e--lower  slot;  b,d,f — upper  slot 


Figure  2  shows  a  non-symmetrical  slotted  wall,  with  differing 
Mach  numbers  and  the  wall  bypass  flow  and  bias  zero  lift  coeffi¬ 
cients.  From  the  figure,  we  can  see  that  when  the  Mach  number  is 
increased,  the  bias  lift  coefficient  produced  on  the  two  wall  sur¬ 
faces  is  also  increased.  Moreover,  the  two  surface  areas  are  larger. 
However,  when  the  bias  quantity  for  the  test  wall  is  produced,  the 
previously  mentioned  interference  becomes  quite  severe. 

Figures  3  and  4  differentiate  the  IAG  and  Xian  wind  tunnel 
upper  and  lower  slotting  pressure  distributions.  In  the  former, 
the  slotting  is  non-symmetrical,  which  causes  the  upper  and  lower 
pressure  distributions  to  be  different.  The  latter  is  symmetrically 
slotted,  and  no  matter  how  large  the  Mach  number  is,  the  upper 
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and  lower  slotting  pressure  distribution  calculation  values  remain 
symmetrical . 


Figure  4.  Northwestern  Polytechnical  University  wind 
tunnel  slotted  wall  pressure  distribution  curve. 

From  the  two  figures,  we  can  see  that  for  a  slotted  wall  area 
aerodynamic  flow  in  experimental  fields  with  slotted  walls  is 
different  from  flow  in  other  domains.  With  different  Mach  numbers 
the  range  of  aerodynamic  flow  in  experimental  stages  will  also 
differ. 
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FINITE  DIFFERENCE  COMPUTATION  OF  THE  FLOW 
AROUND  AIRFOILS  IN  TWO-DIMENSIONAL 
TRANSONIC  SLOTTED  WALL  WIND 
TUNNEL 

Zhang  Xaipir.g 

<  Northwestern  Poiyttchnical  University. 

Abstract 

The  transonic  flow  around  N'ACA  0012  and  RAE  104  airfoils  in  a 
slotted  wail  transonic  wind  tunnel  is  calculated  in  this  paper  with  the 
finite  difference  method. 

A  two-dimensional  small  disturbance  velocity  potential  equation 
is  adopted  in  tnis  computation.  The  transonic  airfoil  wind  tunnels  in 
the  Institute  of  Aerodynamics  and  Gasdynamics  of  the  Stuttgart  Uni¬ 
versity  and  in  the  Institute  ef  Aerodynamics  of  Northwestern  Polytech- 
nical  University  in  Xian  were-  chosen  as  two  computational  examples. 

Only  the  solid  blockage  interference  at  zero  angle  of  attack  is 
calculated  in  this  paper.  The  pressure  distributions  of  the  airfoil  sur¬ 
face  and  the  slotted  wall  along  the  streamwise  direction,  the  addition¬ 
al  lift  coefficient  due  to  the  unsymmetrical  set  up  of  the  model  in 
the  test  section  are  computed. 

The  calculated  results  of  the  N’ACA  0012  and  RAE  104  airfoils  are 
compared  with  the  experimental  results  of  the  Langley  Reseach  Center 
and  those  of  the  National  Physical  Laboratory  in  England  respectively. 
In  general,  the  pressure  distributions  of  the  airfoil  surface  were  sim¬ 
ulated  to  those  of  the  experiments  for  the  selected  Mach  Numbers. 
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